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The  association  of  surfactant  molecules  into  aggregates 
known  as  micelles  gives  them  a  broad  range  of  applications. 
In  spite  of  the  widespread  use  of  this  class  of  chemicals, 
there  is  not  yet  sufficient  scientific  understanding  to 
predict  their  behavior  in  solution. 

A  molecular  thermodynamic  model  has  been  developed  to 
describe  the  formation  of  micelles  in  multicomponent 
surfactant  solutions.   Using  a  hypothetical,  reversible 
seven-step  process,  the  total  free  energy  of  micellization 
is  calculated  by  summing  the  contributions  due  to 
solvophobic  interaction,  mixing,  surface  formation,  confor- 
mational change,  head  group  interactions  and  electrostatics, 
Distributions  of  micelle  sizes  and  compositions  can  then  be 
generated  through  a  set  of  reaction  equilibria.   Where 


possible,  the  free  energy  contributions  are  related  to 
comparable  processes  on  which  experimental  measurements  have 
been  made.   Aggregate  size  distributions  have  been  generated 
from  the  model  for  single-component  solutions  of  nonionic 
surfactants  of  different  chain  lengths  and  at  different 
temperatures  and  solution  concentrations. 

It  has  been  found  that  a  detailed  description  of  micelle 
structure  and  entropy  effects  on  chain  conformation  is 
necessary  to  fully  describe  the  thermodynamics  of  micelle 
formation  without  empirical  parameterization.   To  this  end, 
computer  simulations  of  model  micelles  have  been  conducted 
by  the  molecular  dynamics  method.   Micelles  of  three 
different  head  group  characteristics  and  a  comparable 
hydrocarbon  droplet  have  been  simulated.   A  spherical  shell 
is  used  to  contain  the  aggregate,  providing  estimated 
solvophobic  interactions  with  the  molecules. 

The  simulation  results  reveal  that  internal  structure  of 
the  aggregates  is  relatively  insensitive  to  the  head  group 
characteristics  with  the  greatest  effect  resulting  from  head 
group  size.   While  the  micelles  all  showed  chain  ordering 
and  the  hydrocarbon  droplet  did  not,  the  bond  conformations 
averaged  approximately  71  percent  trans  in  all  cases. 
Results  of  other  simulations  and  experimental  studies  are 


vx 


generally  similar  to  those  of  the  present  work  for  the 
effects  of  chain  length,  aggregate  size,  and  simulation 
technique  on  static  properties. 
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CHAPTER  1 
INTRODUCTION 

There  are  few  classes  of  chemical  species  which  have 
received  more  attention  in  the  scientific  literature,  or 
have  exhibited  a  more  ubiquitous  presence  in  everyday  life, 
than  the  amphiphilic  molecules  known  as  surfactants.   Their 
unique  solution  properties  of  association  and  adsorption  at 
interfaces  give  rise  to  a  broad  range  of  applications.   Long 
the  essential  component  in  detergency  applications,  in  more 
recent  times  surfactants  have  assumed  roles  of  importance  in 
applications  as  diverse  as  enhanced  oil  recovery  and 
phairmacy. 

The  association  of  surfactant  molecules  in  solution  into 
aggregates  known  as  micelles  is  the  primary  attribute  which 
has  garnered  interest  among  the  scientific  community.   Since 
the  discovery  of  micelles  in  solution  (McBain  and  Salmon, 
192  0) ,  many  experimental  and  theoretical  studies  have  been 
conducted,  yet  full  understanding  of  these  systems  of 
molecules  has  not  been  achieved.   True  predictive  capability 
is  not  yet  a  reality. 

This  investigation  takes  a  two-fold  approach  toward  the 
ability  to  make  quantitative  predictions  of  the  behavior  of 
systems  of  surfactants  in  solution.   Since  the  behavior  of 


the  surfactant  solution  is  a  consequence  of  its  thermodynam- 
ics, a  model  set  in  the  framework  of  molecular 
thermodynamics  is  developed  to  describe  the  formation  of 
micelles  in  a  multicomponent  surfactant  solution. 

A  thermodynamic  description  of  a  micellar  system  is 
limited  by  a  precise  knowledge  of  the  structure  of  micelles. 
To  this  end,  computer  simulations  of  model  micelles  by  the 
molecular  dynamics  method  are  conducted.   By  accurately 
modeling  the  forces  present,  pertinent  descriptions  of 
micelle  structure  are  obtained. 

In  Chapter  2  of  this  work,  further  background  on 
micellar  systems  is  given  and  the  development  of  a  model  for 
the  free  energy  change  upon  formation  of  micelles  in  a 
multicomponent  surfactant  solution  is  described.   A  multi- 
step  reversible  process  is  employed  to  generate  contribu- 
tions to  the  total  free  energy  change  due  to  hydrophobic 
interaction,  mixing  effects,  conformational  change,  head 
group  interaction,  and  electrostatic  interaction. 

Some  brief  results  obtained  with  the  thermodynamic  model 
are  presented  in  Chapter  3.   Distributions  of  free  energy 
change  of  micellization  and  aggregate  size  are  calculated 
for  surfactants  of  different  chain  lengths  and  for  solutions 
of  different  temperature.   Due  to  limitations  in  the  scope 


of  this  portion  of  the  investigation,  calculations  were 
carried  out  only  for  cases  of   single-component  solutions  of 
nonionic  surfactants. 

The  computer  simulation  of  surfactant  micelles  is 
described  in  Chapter  4.   Following  a  concise  description  of 
the  molecular  dynamics  method,  its  application  to  the 
simulation  of  micelles  is  discussed.   A  summary  of  the 
computer  simulations  of  the  four  models — three  micelles  and 
one  hydrocarbon  droplet — is  given. 

In  Chapter  5,  the  results  of  analyses  of  the  computer 
simulations  are  presented.   Elements  of  aggregate  internal 
structure,  chain  conformations,  aggregate  shapes,  and 
changes  in  the  aggregate  with  time  are  investigated  for  the 
four  model  aggregates. 

Chapter  6  provides  a  summary  of  the  significant 
conclusions  of  this  work.   Recommendations  are  made  toward 
the  future  progress  of  both  of  the  projects  described  in  the 
previous  chapters. 


CHAPTER  2 
A  MOLECULAR  THERMODYNAMIC  MODEL  OF  MICELLE  FORMATION 

2 . 1  Background 

The  forces  present  in  liquid  solutions  dictate  that 
solution  of  a  polar  solute  in  a  polar  solvent  is  more 
favored  than  is  a  polar  solute  in  a  nonpolar  solvent. 
Similarly,  the  nonpolar  solvent  is  more  accommodating  toward 
a  nonpolar  solute  than  a  polar  one.   Therefore,  molecules 
which  contain  both  polar  and  nonpolar  groups  exhibit  a 
unique  behavior  when  present  in  a  polar  or  nonpolar  solvent. 
Such  molecules,  known  as  surfactants,  will  tend  to  minimize 
the  unfavored  contact  (i.e.,  polar-nonpolar)  while  maximiz- 
ing the  favored  contact  (i.e.,  polar-polar).   At  an 
interface  between  polar  and  nonpolar  liquids,  the 
surfactants  will  penetrate  the  interface  to  achieve  "like" 
interactions  on  both  sides,  reducing  the  "unlike" 
interactions  encountered  in  the  bulk  liquid.   The  polar 
groups  of  a  large  number  of  surfactant  molecules  packed 
closely  at  an  interface  will  repel  each  other,  producing  a 
spreading  pressure  which  reduces  the  interfacial  tension. 

In  the  bulk  liquid,  surfactant  molecules  will  aggregate 
into  structures  known  as  micelles  which  can  afford  much  the 
same  benefits  as  the  interface.   In  the  typical  case,  a 


surfactant  with  a  polar  "head  group"  and  a  nonpolar  "tail," 
when  present  in  sufficient  quantity  in  a  polar  solvent  such 
as  water,  will  form  micelles  having  an  interior  consisting 
of  tails  and  possibly  some  nonpolar  solubilizate  and  a 
surface  comprised  mostly  of  head  groups.   The  head  groups 
remain  in  contact  with  the  water--a  favored  interac- 
tion— while  the  tails  reduce  their  contact  with  the 
water — an  unfavored  interaction.   When  the  solvent  is 
nonpolar,  inverted  micelles  can  form. 

The  ability  of  surfactants  to  adsorb  at  interfaces  and 
aggregate  into  micelles  makes  them  a  very  useful  class  of 
compounds.   The  reduction  of  interfacial  tension  has  many 
applications,  ranging  from  oil  recovery  to  biological 
processes.   In  addition,  micelles  can  solubilize  other 
solutes  in  their  interior,  as  in  drug  delivery  processes, 
and  reactions  can  even  take  place  there,  as  in  emulsion 
polymerization.   The  earliest  and  best  known  use  of 
surfactants,  detergency,  uses  both  aspects  of  their 
behavior.   As  useful  as  these  phenomena  are,  they  are  not 
understood  to  a  degree  that  would  allow  their  full  potential 
to  be  realized.   The  ability  to  predict  behavior  rather  than 
just  explain  it  is  the  goal  of  this  undertaking.   This 
requires  a  knowledge  of  the  thermodynamics  of  surfactant 
phenomena. 


The  formation  of  surfactant  aggregates  in  solution 
instills  a  certain  ambiguity  in  the  description  of  the 
system  by  a  traditional  thermodynamic  formalism.   The 
aggregation  of  surfactant  monomers  into  micelle  structures 
has  been  treated  as  the  formation  of  a  "phase"  (Blankschtein 
et  al . ,  1985;  Kamrath  and  Franses,  1983;  Matijevic  and 
Pethica,  1958)  or  as  a  stepwise  association  "reaction" 
(Tanford,  1974;  Mukerjee,  1972;  Murray  and  Hartley,  1935). 
Although  the  former  description  may  aid  in  visualizing 
certain  aspects  of  micellar  solutions,  the  thermodynamic 
idea  of  a  phase  cannot  be  used  in  a  rigorous  manner.   Its 
requirements  of  continuity  and  homogeneity  are  not  met  by  a 
collection  of  micelles  in  solution  and  a  single  micelle 
cannot  be  treated  as  a  phase  since  its  properties  are 
size-dependent . 

The  treatment  of  micelle  formation  as  reaction 
equilibrium  is  plagued  by  the  lack  of  a  single 
stoichiometry .   Since  a  distribution  of  products  is  formed 
(Void,  1950) ,  one  must  consider  each  micelle  to  be  in 
reaction  equilibria  with  the  dispersed  monomers.   The 
determination  of  the  many  equilibrium  constants  by 
experimental  methods  is  impossible.   Hall  and  Pethica  (1967) 
proposed  using  a  small-systems  thermodynamics  approach  to 


avoid  the  difficulties  of  these  two  treatments.   But  their 
approach  cannot  be  used  with  ionic  systems  and  is  mainly 
formal,  not  lending  itself  to  practical  use. 

The  thermodynamics  of  micelle  formation  remains  an 
interesting  problem.   The  literature  is  abundant  with 
studies.   In  addition  to  the  quantities  of  temperature, 
pressure,  and  composition  which  typically  define  the 
thermodynamic  state  of  a  typical  solution,  the  thermodynamic 
behavior  of  solutions  containing  surfactant  species  can 
depend  on  the  size,  shape,  and  structure  of  the  aggregates 
which  are  formed.   Thermodynamic  properties  have  been 
measured  and  correlated  (Burchfield  and  Woolley,  1984; 
Woolley  and  Burchfield,  1984,  1985).   Aggregate  formation 
has  been  investigated  from  the  points  of  view  of  classical 
thermodynamics  (Moroi  et  al . ,  1984;  Muller,  1973)  and 
statistical  thermodynamics  (Hoeve  and  Benson,  1957;  Owenson 
and  Pratt,  1984) .   Investigations  have  focused  on  size 
distributions  (Ruckenstein  and  Nagarajan,  1975;  Ben-Naim  and 
Stillinger,  1980),  the  role  of  micelle  shape  (Tanford,  1974; 
Israelachvili  et  al.,  1976;  Ljunggren  and  Eriksson,  1984, 
1986;  Eriksson  and  Ljunggren,  1985;  Void,  1985) ,  and  shape 
transitions  (Van  de  Sande  and  Persoons,  1985;  Ikeda,  1984; 
Missel  et  al.,  1983;  Mukerjee,  1977). 
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While  contributing  to  our  understanding  of  the  complex 
nature  of  micelle  formation,  none  of  these  works  produced  a 
practical  model  with  predictive  capabilities.   A  semi-empir- 
ical model  for  the  thermodynamic  properties  of  surfactant 
aggregate  formation  based  on  molecular  thermodynamic 
processes  was  developed  by  Hourani  (1984)  and  was  successful 
at  predicting  thermodynamic  quantities  and  aggregate  size 
distributions  for  systems  of  a  single  surfactant  species  in 
solution.   Benedek  (1985)  developed  a  different  model  in  the 
framework  of  molecular  thermodynamics.   While  it  was 
demonstrated  successfully,  the  use  of  empirical  parameters 
was  more  extensive  than  in  Hourani 's  work.   The  model  of 
Hourani  showed  promise  of  being  extendable  to  multicomponent 
systems  and  of  being  more  closely  related  to  other 
observable  molecular  phenomena.   The  beginnings  of  such  an 
extension  are  given  in  this  chapter. 

2.2  Stoichiometry  and  Reaction  Equilibria 
for  Multicomponent  Micelles 

Since  the  free  energy  of  a  process  is  independent  of  the 
path  chosen  between  the  initial  and  final  states,  Hourani 
proposed  a  process  consisting  of  a  series  of  steps  for  which 
the  free  energy  change  can  be  modeled.  In  this  process,  the 
monomers  were  removed  from  the  solvent  and  placed  in  a 
vacuum  at  their  original  density — a  gaseous  state.  Cavities 
of  excluded  volume  remained  in  the  solvent,  to  be  coalesced 


in  a  subsequent  step.   The  monomer  gas,  considered  ideal, 
was  compressed  to  micellar  density,  and  placed  into  larger 
cavities  which  had  been  formed  in  the  solvent.   Essential  to 
the  modeling  was  the  elimination  and  creation  of  the  solvent 
cavities.   The  counter ions  were  handled  in  the  same  fashion, 
with  the  addition  of  the  necessary  electrostatic  calcula- 
tions.  Extending  Hourani's  molecular  thermodynamic  model  to 
solutions  containing  two  or  more  surfactant  species  required 
the  addition  of  new  steps  and  modification  of  others.   The 
solvent  cavity  steps  were  eliminated  and  the  monomers  are 
removed  to  a  liquid  state  rather  than  the  gaseous.   These 
changes  facilitate  the  handling  of  multiple  components.   The 
development  of  the  multicomponent  model  is  detailed  below. 
Micelle  formation  in  solution  yields  a  distribution  of 
micelle  sizes.   In  addition,  a  multicomponent  surfactant 
solution  has  a  distribution  of  compositions  (Warr  et  al., 
1983;  Scamehorn  et  al.,  1982;  Birdi,  1975;  Moroi  et  al., 
1974,  1975a,  1975b;  Rubingh,  1979;  Clint,  1975).   To 
describe  this,  a  set  of  equilibrium  reactions  can  be  written 
for  the  formation  of  J  micelles  of  distinct  sizes  and/or 
compositions.   For  I  different  surfactant  species,  Z^,  and  K 
counterion  species,  Bj^,  aggregating  to  form  J  micelles,  Z J , 
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(2.1) 


N,jZ,  +  N^jZ^^...^N,jZ,^M,jB,  +  M^jB^^...^M,jB^^Z' 
where  N-j^j  is  the  number  of  monomers  of  species  i  present  in 
the  j'th  micelle  and  Mj-j  is  the  number  of  k  counterions  bound 
to  it.   The  equilibrium  constant  Kj  for  the  formation  of  the 
micelle  is  given  by 

M 

These  are  related  to  the  standard  state  free  energy  of 
micellization  of  the  j'^^  micelle  by 

AG°  =  -RT\nK,  (2.3) 

The  total  number  of  monomers  in  this  micelle  is 

A^;  =  IA^,;  (2.4) 

I 

so  that 

N  ,AQ,°  =  -RT\nK ,  (2.5) 

where  JGj°  is  on  a  per  monomer  basis. 

The  mole  fraction  of  monomers  in  the  j'^h  micelle  is 
found  by  combining  equations  (2.2)  and  (2.5).   For  the 
dilute  concentration  range  of  micellar  solutions,  ideality 


11 
of  the  monomer  solution  can  be  assumed,  and 


^  =  ^"P 


A/,(^-^-I.v,,ln[z,]j-XA^,,ln[B,]-lnC„ 


(2.6) 


where  Cq  is  the  total  solution  concentration 
[  ]  indicates  concentration  of  species 
XJ  is  the  monomer  mole  fraction  in  j'*^'^  micelle 
Nj  is  the  aggregation  number  of  the  j^h  micelle 
Within  the  material  balance  constraint  for  each  species, 

x,  =  Xv,,  (2.7) 

equation  (2.6)  describes  the  distributions  of  micelle  size 

and  composition  in  solution  when  provided  with  the  free 

energy  change  as  a  function  of  the  size  and  composition  of 

the  micelle  formed. 

2 . 3  Estimation  of  Free  Energy  Changes 
The  standard  state  free  energy  of  micellization  is 

calculated  via  the  seven-step  process  shown  in  Figure  2-1. 

The  standard  state  free  energy  of  formation  for  a  single 

micelle  in  solution  is  the  sum  of  the  free  energies  of  the 

seven  steps.   Each  step  is  modeled  as  closely  as  possible 

from  an  observable  phenomenon  of  a  similar  nature. 

Certain  free  energy  terms  are  dependent  on  the  shape  of 

the  micelle.   The  micelle  grows  as  a  roughly  spherical 

aggregate  until  the  additional  volume  of  one  more  monomer 

would  cause  the  radius  of  the  sphere  to  exceed  the  length  of 

the  longest  all-trans  chain.   With  the  addition  of  the  next 
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monomer  the  micelle  must  grow  with  a  nonspherical  geometry 
to  avoid  the  formation  of  a  material-free  core.   Several 
geometries  have  been  proposed:   spherical  dumbells,  oblate 
ellipsoids,  prolate  ellipsoids,  and  spherocylindrical  rods. 
While  Ljunggren  and  Eriksson  (1984,  1986;  Eriksson  and 
Ljunggren,  1985)  have  proposed  that  the  shape  fluctuates 
between  spherical,  rod-shaped,  and  even  disc-shaped.  Void 
(1985)  has  found  little  effect  of  the  particular  geometric 
model  on  the  thermodynamics  of  micelle  formation.   In  this 
work  nonspherical  micelles  are  modeled  as  prolate 
spherocylindrical  rods.   The  derivation  of  micelle  dimen- 
sions and  surface  area  based  on  this  geometric  model  is 
given  in  Appendix  A. 

In  step  1,  a  standard  state  solution  of  surfactants  is 
transformed  into  a  solution  of  hydrocarbon  chains  by 
reversibly  removing  the  head  groups  and  counterions.   Since 
these  are  reversibly  replaced  in  steps  6  and  7,  the  net  free 
energy  change  for  the  mere  removal  and  replacement  of  the 
head  groups  and  counterions  is  zero.   If  no  free  energy 
contributions  due  to  the  replacement  of  head  groups  and 
counterions  are  contained  in  steps  6  and  7,  then 

Z1G,°  =  0  (2.8) 

In  step  2  the  hydrocarbon  solution  is  separated  into  I 
pure  hydrocarbon  liquids  and  the  pure  solvent.   This  is  the 
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reverse  process  of  hydrocarbon  solubility,  so 


=  XI'V,,lnCr  (2.9) 


In  step  3  the  I  pure  hydrocarbons  are  placed  into  J 
ideal  hydrocarbon  mixtures  of  different  compositions  and 
amounts.   For  this  ideal  mixing, 


^7 


IIx,^ln.v,,  (2.10) 


In  step  4  the  J  hydrocarbon  mixtures  are  formed  into  J 
droplets  and  placed  into  the  solvent.   The  free  energy  of 
this  step  is  the  free  energy  of  forming  the  hydrocarbon-sol- 
vent interface.   There  are  both  surface  area  and  curvature 
contributions  to  this  step.   The  expression  for  ^64°,  based 
on  Buff  (1955)  and  Stillinger  (1973) ,  is  the  surface  area  of 
the  droplet,  S,  times  the  planar  interfacial  tension,  y,  of 
the  hydrocarbon  mixture,  corrected  for  curvature: 

-=^  =  — ^  1  -^  (2.1  la) 

RT        NkT\        RJ  ^  ^ 

The  curvature  effect  is  dependent  on  the  parameter,  g.   The 
surface  area  depends  on  the  size  and  shape  of  the  micelle. 
For  the  spherocylindrical  micelle,  a  second  curvature 
parameter,  gc,  is  used  for  the  cylindrical  portion: 


AG^°     2nRy 


RT         NkT 


2R\   1-^Udl  ^' 


R  \         R 


(2.1  lb) 
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This  approximates  the  cylindrical  curvature  effect,  whose 
uncertainty  has  been  discussed  by  Henderson  and  Rowlinson 
(1984) . 

In  step  5  conformational  changes  in  the  hydrocarbon 
chains  are  made  so  that  one  end  of  each  chain  is  at  the 
surface  of  the  droplet.   The  contribution  from  this  step  is 
entirely  entropic  and  may  only  be  estimated.   The  expression 
used  in  this  model  is 

-=^--SJ   ^  (2.12) 

RT  \RJ 

where  S^  is  a  parameter  of  the  model.   The  squared  ratio  of 
chain  length,  1^,  to  micelle  radius,  R,  takes  into  account 
the  very  severe  conformational  restrictions  present  when  the 
micelle  radius  is  much  smaller  than  the  chain  length. 

As  indicated  in  the  discussion  of  step  1,  step  6 
contains  no  contribution  due  to  the  reattachment  of  head 
groups.   The  quantity  JGq     is  the  free  energy  change  due  to 
the  interaction  between  the  head  groups  in  their  positions 
at  the  micelle  surface.   The  head  groups  are  modeled  as 
dipoles.   For  the  ionic  surfactant  species,  a  charged  head 
group  paired  with  a  counter  ion  forms  a  strong  dipole. 
Nonionic  head  groups  exhibit  weak-to-moderate  dipole 
moments.   The  dipole-dipole  interactions  of  adjacent  pairs 
are  summed  for  the  free  energy  contribution  of  this  step. 
The  separation  and  orientation  between  two  dipoles  are 
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dependent  on  the  size  and  shape  of  the  micelle,  with  the 
head  groups  evenly  spaced  over  the  surface  of  the  micelle. 
The  derivation  of  the  head  group  interactions  is  carried  out 
in  Appendix  B.   The  potential  between  a  pair  of  adjacent 
head  groups  is 


b  12 


Dr- 


r    '' 


2R 


(2.13) 


where  n    is  the  dipole  moment 

r  is  the  pair  separation 

R  is  the  micelle  radius 

D  is  the  dielectric  constant  of  the  solvent 
As  indicated  in  the  derivation,  this  form  of  the  potential 
takes  into  account  the  angle  between  adjacent  dipoles  as  a 
function  of  micelle  radius.   The  total  energy  contribution 
from  this  step  is  found  by  summing  the  contributions  from 
the  different  pairs  in  the  manner  described  by  equations 
(B.25)  through  (B.28)  in  Appendix  B. 

Step  7  contains  no  contribution  for  the  replacement  of 
the  counterions  back  into  solution.   The  free  energy  of  this 
step  is  due  to  the  difference  between  the  original  random 
distribution  of  counterions  in  the  micelle-free  solution  and 
the  final  Poisson-Boltzmann  distribution  of  counterions 
around  the  surface  of  the  micelle  with  a  fraction  bound  in  a 
Guoy-Chapman  electrical  double  layer.   The  derivation  of 
Hourani  (1984)  for  the  numerical  solution  of  this  charge 
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distribution  model  is  applicable  here.   The  entire  step  is 
actually  the  process  of  discharging  the  counterions  in  their 
original  distribution,  compressing  them  into  the  bound  layer 
and  final  distribution,  and  then  recharging  the  counterions. 
Therefore  JGy"  contains  an  entropy  contribution  from  the 
compression  and  an  enthalpy  contribution  from  the 
distribution.   The  bound  layer  will  be  populated  with 
dipoles  formed  by  head  group/counterion  bound  pairs  and 
unbound  head  groups.   The  charge  interactions  in  the  bound 
layer  are  included  in  the  dipole  and  point  charge  pairings 
of  the  head  group  term,  JGg". 

2.4  The  Calculational  Technique 
The  free  energy  of  formation  for  a  single  mixed  micelle 
with  nonionic  head  groups  is  obtained  by  summing  the 
contributions  from  steps  1  through  6: 


AG.° 


RT       V^  '  '  '^       NkT 


^"-^J^^M'-^.. 


-SA 


R 


—  1 
RT 


(2.14) 


heads 


To  include  the  presence  of  ionic  head  groups  and 
counterions,  the  contributions  of  step  7  must  be  added  to 
equation  (2.14).   Such  calculations  were  not  carried  out 
here,  so  this  section  will  pertain  only  to  mixtures  of 
nonionic  surfactants.   Table  1  summarizes  the  model's 
variables,  parameters,  and  required  data. 
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Table  1 
Arguments  of  the  Thermodynamic  Model 


Variables 


Data 


Parameters 


T 

C 
C 


tot 
eq 


D 

la 

g 

Sr 


Compositions  of  micelles 

Aggregate  sizes  of  micelles 

Temperature  of  system,  Kelvin 

Composition  of  system 

Number  of  carbons  in  surfactant  chain 

Overall  system  concentration,  moles/liter 

Solubility  concentrations  of 

hydrocarbons,  moles/liter 
Interfacial  tensions  of  hydrocarbon 

mixtures,  dynes/ cm 
Dipole  moments  of  head  groups,  Debyes 
Dielectric  constant  of  water 
Surfactant  chain  length.  Angstroms 
Surfactant  chain  volume.  Angstroms^ 
Spherical  curvature  parameter,  Angstroms 
Cylindrical  curvature  parameter.  Angstroms 
Entropy  of  conformation  parameter 
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The  values  of  chain  length  (Ang.)  and  chain  volume 
(Ang.3)  used  in  the  determination  of  micelle  size  and  shape 
are  calculated  from  Tanford's  correlations  (Tanford,  1972): 

Z,=  1  .265rtc+ 1 -5  (2.15) 

v,  =  26.9nc  +  27  A  (2.16) 

To  facilitate  the  use  of  computer  programs  in  carrying 
out  the  calculations,  correlations  are  used  for  the  required 
physical  data.   The  aqueous  solubility  of  hydrocarbons  used 
in  the  calculation  of  -IG2°  is  obtained  from  a  correlation 
due  to  Leinonen  et  al .  (1971): 


eq 

X     =  exp 


(l-2-X--)^-ln(l-.v- 


(2.17) 


The  parameter  K  is  fit  to  the  solubility  data  of  McAulliffe 
(1966),  Polak  and  Lu  (1973),  and  Sutton  and  Calder  (1974). 
It  is  found  to  be  a  linear  function  of  the  hydrocarbon  chain 
length  for  the  n-alkanes,  but  since  hydrocarbon  solubility 
in  water  exhibits  a  break  at  decane,  two  linear 
relationships  for  K  are  used,  one  for  the  longer  chains  and 
one  for  the  shorter  chains.   Equation  (2.17)  is  solved 
iteratively  for  the  hydrocarbon  solubility,  x^'3. 

For  the  hydrocarbon-water  interfacial  tension  required 
in  the  calculation  of  JG4°,  a  correlation  based  on  the  works 
of  Aveyard  et  al.  (1972),  on  the  surface  tension  of 
hydrocarbons,  and  Jasper  (1972) ,  on  the  surface  tension  of 
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water,  is  used: 

57.868nc+  1  17  .99  -  i  .059nc  +  .  1 768 J7 
y=  1.381 ^ —  (2.18) 

where  the  temperature  is  in  Kelvin  and  the  surface  tension 
is  in  dynes/cm. 

The  value  of  the  dielectric  constant  of  water  used  in 
the  evaluation  of  JG5'  is  given  by 

D  =  252. 422 -.806329  7+  .0007469  7^  (2.19) 

which  is  a  polynomial  fit  of  the  data  of  Owen  et  al.  (1961) 

at  atmospheric  pressure  with  temperature  in  Kelvin.   The 
values  of  dipole  moments  used  in  this  calculation  are 
estimated  as  the  dipole  moments  of  molecules  of  similar 
structure  to  the  head  groups. 

The  three  parameters,  g,  gc,  and  S^;,  are  fit  to  the 
measurable  data  on  the  micellar  system.   These  are  the  mean 
aggregate  size  of  the  micelles  and  the  set  of  mixture 
critical  micelle  concentrations  (CMCs)  at  the  system 
temperature.   Equation  (2.6)  generates  an  I-dimensional 
surface  of  the  aggregate  size  distribution  of  the 
multicomponent  micelles.   This  is  accomplished  by  first 
choosing  values  of  the  temperature,  T,  the  total 
concentration,  C-^ot'  ^'^'^   the  system  composition,  x^.   Then 
for  each  possible  micelle  composition,  x-j^j  ,  equation  (2.6) 
is  evaluated  at  values  of  Nj  from  two  to  infinity  (in 
practice,  until  Xj/Nj  becomes  negligible) .   This  generates 
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the  distribution  of  aggregate  sizes  with  micelle  composi- 
tion.  The  parameters  must  be  chosen  such  that  the 
distribution  meets  the  material  balance  constraints  of 
equation  (2.7) . 

To  find  the  proper  values  of  the  parameters,  the  mean 
size  is  calculated  as 


U 


N  =  ^ (2.20) 


and  the  CMC  is  taken  to  be  the  value  for  which 

lim  ^ =^  =  0.5  (2.21) 

as  put  forth  by  Hall  and  Pethica  (1967) .  Here,  C^  is  the 
free  monomer  concentration  and  Cj^  is  the  concentration  of 
micelles,  calculated  as 

Cm=      °^  (2.22) 

Once  a  set  of  parameters  is  determined  for  a  system,  the 
response  of  the  size  distribution  to  changes  in  any  of  the 
input  variables  can  be  investigated. 

The  program  listing  in  Appendix  C  gives  the  FORTRAN 
source  for  the  interactive  fitting  of  parameters  and 
calculation  of  aggregate  size  distribution  on  systems  of  a 
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single  nonionic  surfactant  species  in  water.   This  is  the 
simplest  application  of  the  model  and  was  used  to  generate 
the  results  presented  in  the  next  chapter. 


CHAPTER  3 
RESULTS  OF  THERMODYNAMIC  MODELING 

3 . 1  Parameter  Estimation 
Three  parameters  of  the  model,  g,  gc,  and  Sq,    must  be 
found  by  fitting  the  mean  aggregate  size  generated  by  the 
model  to  the  experimental  value  at  the  critical  micelle 
concentration  (CMC) .   The  model  output  is  considered  to  be 
at  the  CMC  when  the  total  surfactant  concentration  is  equal 
to  the  experimental  value  of  the  CMC  and  the  derivative  of 
the  monomer  concentration  with  respect  to  the  total 
concentration  satisfies  equation  (2.21),  indicating  that  one 
out  of  two  surfactant  monomers  added  to  the  solution  at  this 
concentration  would  join  a  micelle.   The  model  convergence 
is  quite  sensitive  to  the  values  of  the  parameters,  so  in 
fitting  them  to  the  data,  one  must  either  use  good  initial 
guesses  or  approach  the  values  conservatively.   Failure  to 
do  either  of  these  can  result  in  an  aggregate  size 
distribution  which  has  infinite  values  of  C-j-ot  ^^'^  ^' 
providing  no  useful  information. 

Each  of  the  parameters  has  a  physical  significance  which 
can  aid  in  the  choice  of  the  initial  guesses.   The 
parameters  g  and  gc  are  used  in  the  curvature  corrections  to 
the  planar  interfacial  tension  for  the  spherical  and 
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cylindrical  geometries,  respectively.   Physically,  g  is  the 
thickness  of  the  spherical  interface  in  Angstroms  (Buff, 
1955) .   Although  experimental  values  are  not  available,  it 
is  expected  to  be  positive  and  small  relative  to  the  micelle 
radius.   The  parameter  gc  does  not  have  the  same  physical 
connection  to  the  interface  (Henderson  and  Rowlinson,  1984), 
but  by  comparison  it  is  expected  to  be  positive  and  less 
than  g.   The  parameter  Sq   is  the  (dimensionless)  conforma- 
tional entropy  change  for  a  monomer  joining  a  large  micelle. 
Conformational  contributions  of  the  aggregated  monomers  were 
studied  via  a  statistical  thermodynamic  theory  by  Ben-Shaul 
and  coworkers  (Ben-Shaul,  Szleifer  and  Gelbart,  1985; 
Szleifer,  Ben-Shaul  and  Gelbart,  1985)  and  values  of  the 
conformational  entropy  were  found  to  be  in  the  range  of  -8 
to  -7   for  a  chain  with  seven  bonds.   While  these  values  are 
based  strictly  on  theoretical  considerations  with  no 
experimental  corroboration,  the  parameter  Sq   is  expected  to 
be  of  the  same  sign  and  magnitude. 

In  approaching  the  parameter  fitting  conservatively,  the 
initial  values  of  the  parameters  are  chosen  such  that  the 
aggregate  size  distribution  will  definitely  converge.   That 
is,  the  condition 


'(^ 


lim-^ — ^<0  (3.1) 
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must  be  satisfied.   Equation  (2.6)  defines  the  aggregate 
size  distribution.   In  this  equation,  JGj°  is  the  only  term 
influenced  by  the  parameters.   In  order  for  the  distribution 
to  converge  as  N  becomes  large,  this  free  energy  change  must 
be  greater  than  the  logarithm  of  the  monomer  concentration. 
Overestimating  the  parameters  toward  a  less  negative  free 
energy  change  will  assure  that  the  size  distribution 
converges.   The  parameters  can  then  be  adjusted  toward  a 
more  negative  free  energy  change  as  they  are  fit  to  the 
data. 

The  effect  of  adjusting  the  parameters  can  be  foreseen 
by  analyzing  the  corresponding  terms.   The  parameters  g  and 
gc  affect  the  behavior  of  the  surface  free  energy,  ^64°. 
For  g  <  R  and  gc  <  1q,  the  usual  cases,  increasing  the 
parameter  decreases  the  free  energy,  consistent  with 
equation  (2.11).   In  equation  (2.12)  it  can  be  seen  that 
increasing  the  parameter  S^  decreases  the  conformational 
free  energy,  ^65°. 

The  aggregate  size  distribution  can  be  divided  into  two 
regions.   By  designating  the  value  of  N  for  which  R  =  1q  as 
N-trans/  ^^^   transition  of  the  micelle  geometry  from 
spherical  to  spherocylindrical  defines  the  two  regions  of 
the  size  distribution.   For  N  ^  N-^rans'  '^^®  micelle  is 
spherical  with  radius  R.   For  N  >  N^rans'  "^^e  micelle  is 
spherocylindrical  with  radius  1^  and  length  L.   For  sizes 
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above  N-^rans'  ^^®  conformational  free  energy  is  constant 
with  a  value  of  -Sq.      At  N^rans  "^^^  surface  free  energy 
makes  the  transition  from  being  equal  to  the  spherical 
contribution  to  becoming  dominated  by  the  cylindrical 
contribution  for  large  N.   Figures  3-1  and  3-2  show  the 
effect  of  aggregate  size  and  parameter  values  on  the  surface 
and  conformational  free  energies. 

A  typical  aggregate  size  distribution  generated  for  a 
nonionic  surfactant  is  shown  in  Figure  3-3.   A  peak  occurs 
at  N-trans'  beyond  which  the  fraction  of  aggregate  decreases 
with  increasing  N.   The  total  surfactant  concentration  is 
proportional  to  the  area  under  the  distribution  and  is  thus 
influenced  by  both  the  height  of  the  peak  and  the  slope  of 
the  distribution  above  N-j-rans*   Because  the  distribution 
below  N-trans  rises  so  rapidly,  the  mean  aggregate  size,  n, 
is  dependent  only  on  the  slope  of  the  distribution  above 

^trans • 

The  parameter  Sq,    since  it  contributes  over  the  entire 

range  of  N,  affects  both  C-^ot  ^^^  ^'      "^^^   parameter  gc, 
contributing  only  above  N-j-j-^ns'  ^^^   ^  significant  effect  on 
N  but  only  a  slight  effect  on  C-j-q-j-,  while  g  influences  only 
C-tot  ^^d  ^o't  ^'      This  causes  such  interplay  between  the 
parameters  that  a  range  of  parameter  sets  can  produce 
identical  values  of  C-^ot   ^^^  ^'      "^^^  derivative  constraint. 
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(2.21),  can  also  be  satisfied  by  these  parameter  sets  if 

they  are  fit  at  a  value  of  the  free  monomer  concentration, 

C^,  just  below  the  critical  micelle  concentration. 

The  method  used  to  fit  parameters  consisted  of  first 

choosing  low  values  of  the  parameters,  resulting  in  a  C-f-ot 

essentially  equal  to  C^  and  an  n  insignificantly  larger  than 

N-trans  •   '^^®  parameter  gc  was  increased  until  n  reached  the 

desired  value,  and  then  g  was  adjusted  to  result  in  the 

proper  C-^Q-t-.   This  is  repeated  for  different  values  of  Sq   to 

generate  the  range  of  parameter  sets  fitting  the  data.   The 

experimental  data  used  to  generate  the  parameters  is  given 

in  Table  2. 

3.2  Behavior  of  the  Model  for  Single-Component 

Nonionic  Systems 

Parameter  fitting  was  carried  out  on  systems  of 

different  surfactant  chain  lengths  and  at  different 

temperatures.   For  each  system,  linear  relationships  were 

found  to  exist  between  each  pairing  of  the  three  parameters. 

Corresponding  to  the  manner  in  which  the  fitting  was 

accomplished,  the  curvature  parameters  g  and  gc  can  each  be 

expressed  as  a  straight  line  plotted  against  the  S^ 

abscissa.   Such  a  plot  is  given  in  Figure  3-4.   These 

expressions  are  of  the  form 

g  =  mSc  +  b  (3.2) 

gc  =  nS c-*-  d  (3.3) 
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Table  2 
Data  Used  in  Fitting  Model  Parameters  for  Aqueous  Solutions 

of  Nonionic  Surfactants 

Surfactant 


CioH2l(OC2H4)60H 


Ci2H25(OC2H4)60H 
Ci4H29(OC2H4)60H 


a:  Balmbra  et  al .  (1964) 
b:  Corkill  et  al.  (1961) 
c:  Balmbra  et  al .  (1962) 


T  (Kelvin) 

CMC  (M) 

Mean  Size 

298.15 

9.0E-04 

a 

73   a 

308.15 

6.6E-04 

260 

318.15 

6.4E-04 

640 

298.15 

8.7E-05 

b 

400   c 

298.15 

l.lE-05 

a 

3100   a 

32 


CD 

-i-' 

E 

D 
D 

n 


c 
o 

o 


o 
o 


CM 

(— i 

u 

Q) 

O  SI 

£  -p 
4-J       a; 

^  c  ^^ 

O    S  4J 
Q)    0)    U 


U 


m 
o 


0) 

E 

(3 

O  -H 


O     'C 


0)  ^ 

-P  c  o 

-H 

u-j  >  -n 

O  M  c 

CD  fC 

U  — 

c  in  tri 


c    • 

o   u 

-H  en 

W  CO 

0)  o 

U    rH 

G^   • 

0)    M 
^     I 

CO 
(U  o 

s:  CO 


C    CO    rH 

CU  CM     O 
Q)  -H 

t:  -P  ^H 

iJ   (C    0) 

-p  CD  a 

C   r-i    w 
0)    o 

o  x: 


CO 

>  II 

H    O 

C^  en 


a  c 

•H 


o 

C  O 

o 


J9|9LUDJDd    9Jn:^DAjn3 


I 

n   U 

■H 

^  o 

3  -H 
•H    O    rt! 

fc^  c  > 


<4-l 

o 

0) 
=3 


CM 
I 

ID 

n 

CM 


0)  n 

e    • 

(0 
(0 


I 

II 


33 

The  values  of  the  slopes  and  intercepts  found  for  the 
systems  modeled  are  given  in  Table  3 . 

Many  aspects  of  surfactant  behavior  have  been  found  to 
be  linearly  dependent  on  chain  length  and/or  temperature 
(Rosen,  1978) .   The  information  in  Table  3  indicates  that 
such  relationships  are  also  possible  for  the  interaction  of 
the  parameters  of  this  model.   Though  the  three  chain 
lengths  and  three  temperatures  investigated  cannot  conclu- 
sively indicate  linear  behavior,  they  can  indicate  whether 
this  behavior  is  likely.  Linear  regression  of  the  slopes  and 
intercepts  of  Table  3  versus  chain  length  and  temperature 
resulted  in  the  following  equations  and  correlation 
coefficients: 

m  =  -. 00897+  .8817 

5=  .05497-20.20 

n  =  -. 01377+  1  .421 

d=  .1  1257-45.09 
for  the  10  carbon  nonionic,  and 

m  =  -.1464n,- .3170  /?=.9999  (3.8) 

5=1.418/7,-18.10  /?=.9988  (3.9) 

/7  = -.2201/1,- .4676  R=\.000  (3.10) 

d=  1  .553n,-27.29  R  =  .999\  (3.11) 

at  a  temperature  of  298.15  Kelvin. 


=  .9998 

(3.4) 

=  .9778 

(3.5) 

=  .9999 

(3.6) 

=  .9810 

(3.7) 
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Table  3 

Slopes  and  Intercepts  of  Curvature  Parameter  Dependence 

on  Conformational  Parameter 


"q 

T  (Kelvin) 

g  slope 

g  int. 

gc 

slope 

gc  int. 

10 

318.15 

-1.96163 

-2.80775 

-2 

94330 

-9.43866 

10 

308.15 

-1.86885 

-3.15286 

-2 

80330 

-10.1776 

10 

298.15 

-1.78296 

-3.90518 

-2 

66898 

-11.6881 

12 

298.15 

-2.07048 

-1.32359 

-3 

.10874 

-8.80775 

14 

298.15 

-2.36862 

1.74162 

-3 

.54946 

-5.47762 
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Substituting  (3.4)  through  (3.11)  into  (3.2)  and  (3.3), 

the  temperature  relations  for  CIO  nonionic  are 

g  =  (.0549-  .0089  5c)7-20.20+  .88175^         (3.12) 

gc  =  (.1  125-  .0  1  3725c )T- 45.09+  1  -421  45 ^  (3.13) 

and   the   chain   length   relations   at   298.15   Kelvin   are 

g  =  [  1.418-  .14645c)n,-  1  8  .  1  0  -  .3  1705^  (3.14) 

gc  =  (  1  .553-  .2201  5c  )n^- 27. 29-  .46765  c        (3.15) 
These  relations  are  only  valid  for  values  of  Sq   which,  for 

each  chain  length,  produce  values  of  the  curvature 

parameters  which  are  neither  negative  nor  larger  than  the 

micelle  radius.   In  this  sense,  S^  is  dependent  on  chain 

length,  but  there  remains  a  lack  of  uniqueness  in  the 

parameter  fitting  for  the  single  component  case.   This 

results  from  the  finding  that  the  CMC  and  the  concentration 

derivative  (2.21)  are  not  independent  of  each  other. 

The  free  energy  of  micellization  and  its  various 

contributions,  as  given  by  the  model,  are  shown  in  Figures 

3-5  through  3-8.   These  plots  were  made  for  two  chain 

lengths  at  the  same  temperature,  two  temperatures  with  the 

same  chain  length,  and  two  parameter  sets  at  the  same  chain 

length  and  temperature.   Figures  3-5  and  3-8  show  the  effect 

of  two  different  parameter  sets  for  the  same  system.   While 

the  surface  and  conformational  contributions  are  different 

due  to  the  different  parameters,  the  total  free  energy  of 

micellization  is  the  same  in  the  two  figures.   This  follows 
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from  the  fact  that  the  parameter  sets  were  fit  to  the  same 
concentration  and  mean  size  data  and  generate  identical  size 
distributions . 
3.3  Aggregate  Size  Distribution  and  Concentration  Behavior 
With  the  parameter  sets  determined  by  the  methods 
described  above,  the  aggregate  size  distributions  and  the 
effects  of  surfactant  concentration  were  investigated  for 
the  nonionic  alkyl  surfactants.   The  size  distributions  for 
this  class  of  surfactants  exhibit  a  rapid  increase  in  the 
spherical  region  of  N,  peaking  at  N-t-rans*  ^^   ^^^ 
nonspherical  region,  the  distribution  decreases  very  slowly 
through  large  values  of  N,  producing  the  large  mean 
aggregate  sizes  which  have  been  measured  for  the  nonionic 
surfactants . 

The  size  derivative  of  the  distribution  in  the 
nonspherical  region  is  negative.   With  increasing  free 
surfactant  monomer  concentration,  it  becomes  less  negative, 
resulting  in  a  higher  mean  aggregate  size  and  total 
surfactant  concentration.   The  effect  of  surfactant  concen- 
tration on  the  modeled  size  distribution  and  mean  aggregate 
size  of  the  C12  nonionic  surfactant  is  shown  in  Figures  3-9 
and  3-10.   Below  the  CMC,  the  model  calculates  a  mean 
aggregate  size  approximately  equal  to  N^rans'  present  as  the 
initial  plateau  in  figure  3-10.   This  value  is  meaningless, 
however,  since  the  concentration  of  micelles  in  this  range 
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of  surfactant  concentration  is  insignificant.   The  result  is 
merely  the  ratio  of  two  extremely  small  numbers,  since  the 
model  always  produces  a  finite  value  for  the  concentration 
of  micelles  (equation  (2.6)). 

Measurement  of  properties  such  as  surface  tension,  which 
depend  on  the  free  monomer  concentration,  show  that  beyond 
the  CMC  the  free  monomer  concentration  is  virtually 
unaffected  by  an  increase  in  the  total  surfactant 
concentration.   The  relationship  between  the  free  and  total 
surfactant  concentrations,  as  predicted  by  the  model,  is 
exhibited  in  Figure  3-11.   At  total  concentrations  lower 
than  the  CMC,  the  majority  of  any  surfactant  added  to  the 
solution  exists  as  free  monomers,  while  at  concentrations 
higher  than  the  CMC  the  bulk  of  any  additional  surfactant 
micellizes.   The  definition  of  the  CMC  used  in  this  work,  as 
expressed  in  equation  (2.21),  quantifies  the  break  between 
these  two  behaviors  of  the  solution.   Figure  3-12  shows  the 
concentration  derivative  of  equation  (2.21)  for  the  data  of 
Figure  3-11  with  the  CMC  identified. 

3.4  Chain  Length  and  Temperature  Effects 

It  is  known  that  the  hydrophobic  nature  of  the  alkyl 
chain  increases  with  its  chain  length,  as  evidenced  by  the 
decreasing  solubility  in  water  of  the  increasingly  long 
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chains  (McAulliffe,  1966).   With  alkyl  surfactants,  this 
increasing  hydrophobicity  promotes  micellization.   While  no 
significant  aggregate  formation  is  observed  for  chains  of 
six  carbons  or  less,  increasing  chain  length  brings  about 
the  formation  of  larger  micelles  at  lower  critical  micelle 
concentrations  (Table  2). 

The  elements  of  the  model  which  depend  on  chain  length 
are  the  hydrophobic  free  energy,  taken  directly  from  the 
aqueous  solubility,  the  curvature  parameters,  accounting  for 
the  difference  in  micelle  radius,  and  the  conformational 
free  energy.   Equations  (3.14)  and  (3.15)  show  the 
relationship  of  the  parameters  to  the  chain  length.   The 
aqueous  solubility  of  the  alkane,  which  is  chain-length, 
dependent,  is  used  by  the  model  in  the  fitting  of  the 
parameters.   The  dependence  of  the  aggregate  size  distribu- 
tion and  free  energy  change  on  the  surfactant  chain  length 
is  shown  in  Figures  3-13  and  3-14.   The  size  distribution 
broadens  with  increasing  chain  length,  yielding  a  much 
larger  mean  aggregate  size  with  a  less  significant  increase 
in  the  CMC.   The  free  energy  of  micellization  becomes  more 
negative  with  increasing  chain  length,  showing  that  micelle 
formation  by  surfactants  with  more  hydrophobic  area  is  more 
thermodynamically  favored. 
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Just  as  the  aqueous  solubility  of  an  alkane  decreases 
with  increasing  temperature,  the  CMC  of  the  corresponding 
surfactant  decreases  and  the  mean  aggregate  size  increases. 
Again,  by  incorporating  the  solubility  data  into  the  model, 
the  proper  temperature  dependence  is  given  to  the 
parameters.   The  temperature  dependence  of  the  parameters  is 
in  the  form  of  equations  (3.12)  and  (3.13).   The  dependence 
of  the  aggregate  size  distribution  and  free  energy  change  on 
temperature  is  shown  in  Figures  3-15  and  3-16.   The  effect 
of  increased  temperature  is  similar  to  that  of  increasing 
the  surfactant  chain  length:   a  broadened  size  distribution 
and  a  more  negative  free  energy  of  micellization.   But  while 
the  first  increase  of  ten  degrees  Kelvin  has  a  much  more 
pronounced  effect  on  the  free  energy  change  (and  hence  the 
size  distribution)  than  the  next  increase  of  ten  degrees, 
the  mean  aggregate  size  is  more  affected  by  the  latter. 
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CHAPTER  4 
MOLECULAR  DYNAMICS  SIMULATION  OF  SURFACTANT  MICELLES 

4 . 1  Background 
The  development  of  a  predictive  nodel  of  micelle 
behavior,  such  as  that  of  Chapter  2,  requires  a  knowledge  of 
the  structure  of  the  aggregate  in  order  to  describe  its 
thermodynamic  properties.   The  calculation  of  a  surface  free 
energy  contribution  necessitates  certain  assumptions  about 
micelle  size,  shape,  and  surface  roughness.   The  solvent 
interface  must  be  modeled  and  solvent  penetration  must  be 
estimated.   Calculation  of  the  conformational  contribution 
to  the  free  energy  calls  for  details  of  the  molecular 
conformations  in  the  free  and  aggregated  states.   The  head 
group  contribution  depends  on  the  average  positions  of  the 
head  groups . 

^Tiile  experimental  results  are  the  most  desirable  basis 
for  thermodynamic  models,  measurements  of  the  structure  of 
micelles  are  limited  in  the  information  they  can  provide. 
Therefore,  current  models  of  micelle  behavior,  including 
that  derived  in  Chapter  2 ,  are  based  primarily  on 
theoretical  considerations. 

The  most  informative  techniques  for  studying  micelles  in 
solution  are  spectroscopic.   Light  scattering,  Nf-ER,  and 
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small-angle  neutron  scattering  (SANS)  have  been  used 
extensively  to  study  micelles.   However,  results  of  these 
experiments  often  conflict,  as  in  the  question  of  whether 
there  is  solvent  penetration  into  the  micelle  core,  and 
their  description  of  the  micelle  surface  and  internal 
structure  is  incomplete  (Tabony,  1984). 

Quasi-elastic  light  scattering  has  been  used  effectively 
to  determine  micelle  size  and  shape  (Corti  and  Degiorgio, 
1981a,  1981b)  and  Raman  light  scattering  has  been  used  to 
investigate  chain  conformations  (Kalyanasundaram  and  Thomas, 
1976) .   The  absence  of  solvent  from  the  micelle  core  was 
determined  with  NMR  techniques  (Mitra  et  al . ,  1984),  though 
the  opposite  had  been  reported  earlier  (Menger,  1979) .   The 
most  promising  experimental  means  available  of  studying 
micelle  structure  is  SANS.   It  has  been  used  to  study  the 
micelle  surface  (Hayter  and  Zemb,  1982) ,  solvent  penetration 
(Tabony,  1984;  Cabane  et  al.,  1985),  and  chain  conformations 
in  the  micelle  interior  (Bendedouch  et  al . ,  1983a,  1983b). 

The  spectra  generated  by  the  above  methods  must  be 
carefully  interpreted  to  yield  correct  information  about  the 
micelle  surface  and  interior.   It  is  in  this  interpretation 
of  the  results,  particularly  in  accounting  for  intermicellar 
effects,  that  inconsistencies  arise  (Cabane  et  al.,  1985). 
Computer  simulation  is  an  alternative  means  to  answer 
questions  raised  by  the  experimental  measurements.   For  a 
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given  model  of  surfactant  and  micelle  behavior,  simulation 
can  provide  exact  quantitative  descriptions  of  micelle 
structure. 

Simulations  of  micelles  have  been  conducted  by  both 
Monte  Carlo  (Haan  and  Pratt,  1981a,  1981b;  Owenson  and 
Pratt,  1984)  and  molecular  dynamics  (Haile  and  O'Connell, 
1984;  Woods  et  al.,  1986;  Jonsson  et  al.,  1986)  methods. 
Thus  far,  results  of  the  simulation  studies  have  been 
promising.   Their  contribution  to  a  better  understanding  of 
micelle  structure  has  been  limited  only  by  the  uncertainties 
about  the  model  micelle. 

Monte  Carlo  calculations  are  carried  out  on  a  lattice  of 
sites.   Movement  of  a  molecular  segment  from  one  site  to 
another  is  allowed  or  disallowed  based  on  the  energy  change 
associated  with  the  move.   The  restriction  of  segment 
positions  to  the  lattice  sites  limits  the  possible 
conformations  available  to  the  molecules.   To  allow  more 
possible  conformations,  the  number  of  lattice  points,  and 
hence  the  computing  effort,  must  be  increased.   This 
restricts  how  far  such  a  simulation  can  go  in  accurately 
representing  reality;  in  addition,  the  Monte  Carlo  method 
precludes  the  opportunity  to  study  the  dynamics  of  a  system 
by  simulation. 
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The  molecular  dynamics  method  provides  a  continuum  of 

sites  for  molecular  positions  and  is  limited  only  by  the 

accuracy  of  the  forces  built  into  the  model.   In  the  present 

work,  molecular  dynamics  models  were  developed  based  on  the 

work  of  Haile  and  coworkers  (Haile  and  O'Connell,  1984; 

Woods  et  al . ,  1986)  to  investigate  the  effects  of  surfactant 

chain  length,  head  group  mass,  and  head  group  size. 

4.2  The  Molecular  Dynamics  Method 
The  molecular  dynamics  approach  to  computer  simulation 

consists  of  summing  the  forces  on  each  particle  in  the 

system  and  solving  Newton's  equations  of  motion  for  the 

resultant  positions  and  momenta  (Haile,  1980) .   The  problem 

of  constructing  a  realistic  model  system  then  becomes  one  of 

supplying  appropriate  models  for  the  intermolecular  and 

intramolecular  potentials.   It  is  assumed  that  a  given 

particle's  interactions  with  its  neighbors  are  pairwise 

additive.   Thus,  the  potential  between  two  particles  at  time 

t  is  of  the  form 


/:/(r,^(0)  =  t^ir,(0.r/Oj  (4.1) 
The  force  between  the  two  particles  is  given  by 

?(>„(0]  =  -V^(r,/0]  (4.2) 

and  the  total  force  on  the  particle  of  interest  is  the  sum 

of  the  forces  resulting  from  interactions  with  all  other 

particles: 

?,(0  =  I?(r„(0)  (4.3) 
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The  resulting  acceleration  and  velocity  of  the  particle  are 
obtained  from  Newton's  Second  Law  of  Motion,  relating  force, 
mass,  and  acceleration: 

PSn  =  m,^-^  (4.4) 

The  set  of  second-order  differential  equations  for  all 
of  the  particles  in  the  system  is  solved  at  each  time  step 
of  the  simulation.   The  means  employed  in  this  work  for 
determining  the  new  positions  of  all  of  the  particles  from 
the  forces  at  the  previous  time  step  is  a  fifth-order 
predictor-corrector  algorithm  (Gear,  1971) .   This  algorithm 
consists  of  three  procedures  which  are  repeated  at  each  time 
step.   The  position  and  its  first  five  time  derivatives  at 
time  (t+Jt)  are  predicted  for  all  of  the  particles  by 
Taylor's  series  expansion  of  their  values  at  time  t: 

dr,(0     dV,(0(JO'  d\,(t)(iAt)^ 
'^    dt  dt^         2!     dt^         3! 

^  d'r,(t)(Aty  ^  dV.(0(zlO^ 
dt'        4!  ^  dt^        5' 

—  rAt  +  At)= ^^  + ^(JO  + ^- — 

dt    '^     ^    dt  dt^      ^      ^  dt^         2! 

dV,(0(z]0^  dV,(0(JO'* 
dt'        3^^~d?    4V~ 


(4.5) 


(4.6) 


dt'''^    ^       ^         dt^      "   df^   ^  ^"  dt'         2!   ^  dt^         3! 
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d«^  '      ^    dt^  dt'     ^      ^        dt^        2! 


3r,(f^zJ0=-^:T-"— nr-C^O^   ,,s  \.:  (4-8) 


d^         dV,(0  dV,(0 

— -r,(t  +  At)  = -— + 

dt^  dt"  dt 


rXt  +  At)= -— + ^-(Z!0  (4.9) 


d^  cfV,(0 

--r,(f  +  ZlO= i—  (4.10) 

dt^  dt^ 

The  forces  Fi(t+zit)  are  then  evaluated  at  the  predicted 
positions.   Finally,  the  predicted  values  of  the  position 
and  its  derivatives  are  corrected  according  to  the  error 
between  the  acceleration  predicted  by  the  series  expansion 
and  the  acceleration  calculated  from  Fj^(t+^t).   These 
calculations  are  carried  out  by  the  subroutines  PREDCT, 
EVAL,  and  CORR  in  Appendix  D. 

The  simulations  described  in  this  work  were  carried  out 
at  constant  temperature.   The  temperature  of  the  system  is 
related  to  the  velocities  of  the  particles  by 

7  =  — *— Vm.i/f  (4.11) 

The  velocities  of  all  particles  are  scaled  at  each  time  step 
so  that  the  temperature  is  maintained  at  a  constant  value. 
Subroutine  EQBRAT  in  Appendix  D  accomplishes  this  task. 

To  start  a  simulation,  values  of  the  positions  are 
required  to  evaluate  the  initial  forces.   Initial  values  of 
the  five  derivatives  of  position  are  required  by  the  fifth 
order  predictor-corrector  algorithm.   In  these  simulations, 
initial  positions  were  assigned  according  to  a  lattice  of 
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sites  and  random  velocities  were  assigned  the  particles 
(subroutines  INTPOS  and  INTVEL) .   The  momenta  were  scaled  so 
that  there  was  no  net  linear  momentum  of  the  system.   The 
initial  accelerations  are  calculated  by  equation  (4.4)  and 
the  third  and  higher  derivatives  of  position  are  assigned 
initial  values  of  zero,  since  there  is  no  way  of  evaluating 
them.   The  algorithm  recovers  rapidly  from  this  initial 
condition,  arriving  at  proper  values  of  all  derivatives 
within  10-20  time  steps. 

Since  intermolecular  potentials  are  significant  over  a 
pair  separation  which  is  usually  small  relative  to  the 
dimensions  of  the  molecular  system,  computing  effort  can  be 
reduced  by  employing  truncated  potentials  (Haile,  1980)  and 
maintaining  a  neighbor  list  (Verlet,  1967) .   In  a  simulation 
where  the  majority  of  possible  pair  separations  are  so  large 
that  the  corresponding  pair  potentials  will  be  very  small, 
it  is  more  efficient  not  to  calculate  the  potential  for  pair 
separations  which  are  beyond  a  certain  value  and  do  not 
contribute  significantly  to  the  system  properties.   Instead, 
these  contributions  are  neglected  during  the  simulation  and 
a  correction  is  made  to  the  system  properties.   A  cutoff 
distance,  r,-.,  is  set  such  that  the  potential  can  be 
neglected  for  r-^j  >  r^.   The  intermolecular  forces  are 
truncated  at  r^  and  shifted  vertically,  going  to  zero 
smoothly  at  r^   (Nicolas  et  al.,  1979).   The  truncated 
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potential  is  obtained  from  this  shifted  force  by  equation 
(4.2).   This  eliminates  a  step  change  in  the  potential  and 
force  at  r-j^j  =  r,-.,  which  could  have  a  disruptive  effect  on 
the  energy  conservation  of  the  simulated  system.   The 
general  equations  of  a  shifted  force,  Fs(r),  and  its 
potential,  Us(r),  are  given  below: 


f's(r)  = +  r 

dr 


^s(0  =  ^(/-)-f/(rJ-(r-r, 


dr 

dUjr) 
dr 


r<r,  (4.12) 


r<r^  (4.13) 


By  checking  the  value  of  r^j  and  bypassing  the  calculation 
of  F(rij\  for  r^j  >  r^,  a  significant  savings  in  computing 
effort  can  be  realized. 

Further  savings  can  be  made  by  not  checking  pairs  which 
have  a  low  probability  of  r-^j  <  r^.   This  is  accomplished 
through  the  use  of  a  neighbor  list.   A  particle  (i)  in  the 
system  is  surrounded  by  a  sphere  of  radius  r^.   Other 
particles  (j)  inside  this  sphere  will  interact  with  it, 
since  r^j  <  r^.      From  one  time  step  to  the  next,  particles 
will  enter  and  exit  this  sphere.   Over  some  short  period  of 
time,  all  of  the  particles  which  have  interacted  with  the 
central  particle  will  lie  inside  a  larger  sphere  whose 
radius  is  r^isf  ^^ '    over  this  brief  period  of  time,  only 
those  particles  inside  the  larger  sphere  were  checked,  all 
of  the  particles  interacting  with  the  central  particle  would 
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be  found  without  checking  all  of  the  other  particles  in  the 
system.   Through  a  neighbor  list  a  record  is  maintained  of 
all  the  "neighbors"  within  a  distance  r^ist.  of  each  particle 
in  the  system.   Only  these  neighbors  are  checked  for  r3_j  < 
r^.   The  neighbor  list  is  updated  periodically  based  on  the 
value  of  rj^j^g-t^  and  the  dynamics  of  the  system. 

The  molecular  dynamics  method  represents  a  massive 
computing  effort  for  any  system  large  enough  to  be  of 
interest.   However,  using  the  techniques  described  above  to 
improve  efficiency,  and  with  the  development  of  supercomput- 
ers and  parallel  processors,  it  has  become  a  practical  tool 
for  quantitatively  studying  molecular  behavior. 
4.3  The  Model  Surfactant  Molecule 

The  surfactant  molecule  studied  in  this  work  has  a 
linear  alkyl  chain  of  eight  carbon  groups.   All  of  the 
groups  are  considered  to  be  methylenes.   A  polar  head  group 
is  attached  to  one  chain  end;  head  groups  of  different  sizes 
and  masses  were  used  in  the  simulations.   The  bond  lengths 
and  angles  along  the  chain  are  those  of  a  normal  alkane. 
These  are  given  in  Figure  4-1. 

Three  types  of  intramolecular  interaction  take  place  in 
the  model  surfactant  molecule: 

1.  Bond  vibration 

2 .  Bond  bending 

3 .  Bond  rotation 
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1.539  A 


Head 


Group 


Figure  4-1.   The  model  surfactant  molecule,  consisting  of  an 
eight-member  alkane  chain  attached  to  a  polar  head  group, 
shown  in  the  all-trans  conformation.   Bond  angles  and 
lengths  of  the  alkane  chain  are  indicated.   The  bond  to  the 
head  group  is  head  group  dependent  in  the  different 
simulations  conducted. 
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The   bcr.d   vibra'cicr.  pcwsr.-ial    used    is   zr.a.z   cf   vreber    (1978), 
taken   fron  a   simulation    cf   n-butane.       It    is    a   harronic 
potential   aljcut   z'~e   equilibrium  bond   length.      Tr.e  -czzer.zia'. 
eind   the    forces    or.   ~he   rvo    groups    are   giver,   by 


U,{b,]  =  ^y,{b,-bS  (4.14) 


F'Abi]  =  -yAb,-b,]^  (4.15) 

bi 

FlAb^-yAb^-t'o]^^  (-H6) 

where  b^  is  the  bond  lerrth  between  groups  i  and  i-1,  b^  is 
the  eq-uilibriun  bond  lengrh,  and  ■• ,  is  the  force  constant. 

The  bond  bending  pciential  is  also  taken  fro-  the 
simulation  cf  Weber  (1978).   It  is  a  hamonio  potential 
about  the  cosine  cf  ~he  equilibriim  bond  angle  wLz'r.   forces 
on  the  three  groups  foming  -he  bond  angle.   For  the  angle  5 
f erred  by  groups  i,  i-1,  and  i-2, 

lfAo,)  =  ^yJcose,-cose,]^  (4.17) 

F.{e,)'-yJcosG,-cose,  -\^^-^^-cose,\  (4.18) 

-b    f     \         f  [  I   (b.^  1  b,     "^ 

1    (b,  b,*i       w 

—^— COS0,  (4.19) 

b.-l  \b,    b.-l     7j  ^ 
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?!'.2(0j  =  y6(cos0„-cos0,)^f^  +  ^cos0,J  (4.20) 

The  model  surfactant  molecule  is  comprised  of  eight 
bonds  connecting  its  nine  groups.   Rotation  about  any  of  the 
inner  six  bonds  results  in  a  change  in  the  conformation  of 
the  molecule  and  a  change  in  its  potential  energy.   A  set  of 
four  groups  in  sequence  along  the  chain  molecule  contains  a 
dihedral  angle  formed  by  the  three  bonds  connecting  the 
groups.   This  angle  is  a  measure  of  rotation  about  the 
middle  bond,  and  the  potential  energy  is  a  function  of  the 
dihedral  angle.   The  bond  rotation  potential  chosen  for 
these  simulations  is  that  of  Ryckaert  and  Bellemans  (1975) : 

U(<p)  =  y,(\  A\6-  I  .462  cos  0  -  1  .578  cos^  0  +  0  .368cos^0 

+  3.156cos'*0  +  3.788cos^0)         (4.21) 
where  0  is  the  dihedral  angle  in  radians.   The  bond  rotation 

potential  is  illustrated  in  Figure  4-2.   The  force  on  each 

of  the  four  groups  resulting  from  this  potential  and 

equation  (4.2)  has  been  derived  by  Woods  (1985)  and  the 

resulting  calculations  are  carried  out  by  subroutine  RYFOR 

in  Appendix  D. 

The  intramolecular  potentials  account  for  all  of  the 

interactions  among  adjacent  groups  on  a  molecule  as  they  are 

taken  two  (vibration) ,  three  (bending) ,  or  four  (rotation) 

at  a  time.   Two  groups  on  the  same  molecule  which  are 

separated  by  more  than  three  intervening  groups  do  not 
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interact  through  any  of  these  three  potentials,  but  may 

experience  the  intermolecular  potential  if  they  become 

sufficiently  close  to  one  another. 

4.4  The  Model  Micelle 
The  micelle  used  in  the  simulations  contains  24 

surfactant  molecules.   The  individual  groups  of  the 

molecules  interact  through  intermolecular  potentials,  and 

the  micelle  is  surrounded  by  a  spherical  shell  which  models 

the  solvophilic  effect  on  the  head  groups  and  the 

solvophobic  effect  on  the  chains.   The  intermolecular 

interactions  which  take  place  in  the  micelle  are 

1.  Head  group-shell 

2.  Alkane  group-shell 

3 .  Head  group-head  group 

4 .  Alkane  group-alkane  group 

5.  Head  group-alkane  group 

The  shell  does  not  attempt  to  explicitly  model  the 
solvent  molecules.   To  do  so  greatly  increases  the  computing 
requirements  of  the  simulation  and  a  previous  attempt  at 
such  explicit  modeling  (Jonsson  et  al .  ,  1986)  did  not  show 
it  to  be  a  clear  improvement  over  the  present  poten- 
tial-field model.   Rather,  the  solvophilic  and  solvophobic 
natures  of  the  head  group  and  alkane  group  are  modeled  in  a 
simple  manner  intended  to  keep  the  head  groups  close  to  the 
surface  of  the  micelle  and  prevent  the  chains  from  extending 
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out  from  the  micelle.   This  simulates  the  minimum 
chain-solvent  contact  of  accepted  experimental  results 
(Tabony,  1984) . 

The  head  group  interaction  with  the  shell  is  through  a 
harmonic  potential  about  an  equilibrium  radial  position: 

u[r,]  =  y,[r,-rS  (4.22) 

?('-,)  =  -2y,(r,-rJ^  (4.23) 

where  r-[  is  the  position  of  the  head  group  and  rQ  is  the 
equilibrium  position.   The  force  constant  is  y^. 

The  alkane  group  interaction  with  the  shell  is  through  a 
repulsive  potential: 

U[n.)  =  4'f-T  (4.24) 

?(-.J=12e[-^^'-  (4.25) 

where  r^^  is  the  quantity  r^-r^,  r^  being  the  radial 
position  of  the  spherical  repulsive  shell,  and  r^j  and  e  are 
the  radius  and  energy  of  minimum  potential  for  the 
alkane-alkane  intermolecular  potential. 

The  interaction  between  two  head  groups  is  described  by 
a  potential  comprised  of  both  dipole  (r~^)  and  soft-sphere 
(r~12)  repulsions: 


^('-.J 


3    /  ^    \\2 
hh    \  I    '   hh 


r^J         V^.y 


(4.26) 


67 


This  potential  is  designed  to  spread  the  head  groups  apart 
on  the  exterior  of  the  micelle.   The  force  between  a  pair  of 
head  groups  is  obtained  by  applying  equation  (4.2)  to 
equation  (4 . 26) . 


?(rj  = 


3rL  12r'' 


hh 


13 


I2I 

r,, 


(4.27) 


Applying  equations  (4.12)  and  (4.13)  to  use  the  shift- 
ed-force  potential  for  rj^j  <r(-., 


Us[r,, 


hh 


^  hh 


12 


12 


-4 


-13 


^  hh 


13 


(4.28) 


^s(^.;) 
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hh 
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(4.29) 


The  intermolecular  potential  between  two  alkane  groups 
is  a  (6,9)  form  of  the  Mie  (m,n)  potential  (Reed  and 
Gubbins,  1973) : 


^(r,J  =  e 


fir,J=186 


^Xr.y-M'-A  -^['- 


20 


-  li 


'    m        '    m 


+  21 


I1Z 


/^ 


10 


(4.30) 


(4.31) 


(4.32) 
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^■(-'=-{^^^-^[(^)'-(7 


■0   ,  J. 

'—  (4.33) 

This  potential  is  also  used  for  the  interaction  between  a 
head  group  and  an  alkane  group,  except  that  the  radius  of 
minimum  potential,  r-^,    is  adjusted  to  account  for  the 
difference  between  the  diameter  of  the  head  group  and  that 
of  the  alkane  group: 

^h^ead-a.^ane^^-^^^^^J  (4.34) 

Figure  4-3  illustrates  the  intermolecular  potentials  used  in 
the  model  micelle,  and  Table  4  gives  the  values  of  the 
parameters  used  in  the  intramolecular  and  intermolecular 
potentials. 

The  model  micelle  is  assembled  by  initially  placing  the 
24  molecules  with  their  head  groups  evenly  spaced  over  the 
surface  of  a  sphere  of  twice  the  expected  micelle  radius. 
Each  chain  is  directed  radially  inward,  in  the  all-trans 
conformation.   With  the  bond  rotation  force  constant  reduced 
to  one-tenth  of  its  normal  value  to  facilitate  chain  packing 
during  startup,  the  simulation  is  begun.   The  radius  of  the 
confining  sphere  is  reduced  every  ten  time  steps  until  the 
appropriate  radius  is  achieved.   This  is  the  value  that 
would  give  the  density  of  the  analogous  liquid  alkane, 
adjusted  to  a  pressure  which  fluctuates  about  zero  with  an 
amplitude  of  less  than  ten  atmospheres.   After  this,  the 
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Figure  4-3.   The  intermolecular  potentials  for  segment-seg- 
ment (6,9),  segment-shell  (12),  and  head-head  (3,12).   A 
value  of  419  Joule/mole  was  used  for  e;  both  r-^   and  T-^ii   were 
assigned  values  of  4  Angstroms. 
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Table  4 
Parameter  Values  for  Potentials  Used  in  Simulations 
Potential        Parameter     Value  Units 


Bond  vibration  ^ 

y« 

9 

25xl05 

a 

Joule/A^/mole 

bo 

1.539 

a 

Angstrom 

Bond  bending  + 

yt 

1.3xl05 

a 

Joule/mole 

do 

112.15 

a 

degree 

Bond  rotation  "•" 

Vr 

8313 

b 

Joule/mole 

Head-shell 

Vh 

785 

c 

Joule/A^/mole 

ro 

4 

c 

Angstrom 

Segment-shell 

e 

419 

Joule/mole 

rm 

4 

Angstrom 

Head-head 

e 

419 

Joule/mole 

r^h 

* 

Angstrom 

Segment-segment 

e 

419 

a 

Joule/mole 

r,n 

4 

a 

Angstrom 

Head-segment 

e 

419 

Joule/mole 

h-a 

^( 

r^+r^J 

Angstrom 

a  Weber  (1978) 

^   Ryckaert  and  Bellemans  (197  5) 

c  Woods  (1985) 

*  Parameter  varies  with  head  group  size. 

+  See  also  Table  6. 
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rotational  force  constant  is  restored  to  its  desired  value 
and  the  system  is  allowed  to  equilibrate.   When  there  is  no 
net  drift  in  the  energy  of  the  system  or  in  the  fraction  of 
trans  bonds,  equilibrium  is  considered  to  have  been  reached. 
The  simulation  is  continued,  periodically  saving  the 
positions  and  velocities  of  all  of  the  groups  for  subsequent 
analysis. 

The  time  step  used  in  simulation  must  be  small  enough  to 
track  the  motions  of  the  molecules  and  maintain  stability, 
yet  large  enough  to  avoid  unnecessary  computing  effort.   The 
ideal  time  step  can  only  be  found  by  trial,  which  was 
conducted  on  a  simulated  liquid  alkane  system  by  Weber 
(1978).   That  result,  Jt  =  .002  psec,  is  the  basis  for  the 
time  steps  used  in  the  present  work. 

4.5  Summary  of  Computer  Simulations 

The  following  simulations  were  carried  out  on  the  model 
micelle: 

1.  Head  group  mass  and  size  same  as  chain  segment. 

2.  Head  group  mass  greater  than  segment,  size  the  same. 

3.  Head  group  mass  and  size  greater  than  segment. 
In  addition,  a  fourth  simulation  was  carried  out  on  a 
hydrocarbon  droplet  of  24  nine-segment  molecules  by  setting 
the  head  group  to  the  size  and  mass  of  a  methylene  and 
replacing  head-head  and  head-shell  interactions  with  segment 
interactions.   From  these  simulations  the  independent 
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effects  of  head  group  size,  mass,  and  solvophilic  nature  can 
be  observed,  and  by  comparison  with  comparable  simulations 
of  dodecyl  surfactant  micelles  (Woods  et  al.,  1986)  the 
effect  of  chain  length  can  be  evaluated.   A  brief  summary  of 
the  simulations  is  given  in  Table  5. 

To  carry  out  the  first  two  simulations,  24  model 
surfactant  molecules  were  placed  with  their  head  groups  on  a 
sphere  of  radius  24  Angstroms  and  their  chains  directed 
radially  inward  in  the  all-trans  configuration.   The  groups 
were  given  initial  velocities,  the  rotational  force  constant 
was  reduced  to  one-tenth  of  its  normal  value,  and  the 
simulation  was  begun.   Every  ten  time  steps,  the  radius  of 
the  confining  sphere  was  decreased  by  .01  Angstrom.   This 
scaling  down  of  the  model  micelle  was  continued  until  the 
radius  of  the  spherical  shell  reached  a  value  of  12.00 
Angstroms,  slightly  higher  than  the  radius  of  an  analogous 
hydrocarbon  droplet.   Ten-thousand  time  steps  were  run  at 
this  radius  to  allow  the  micelle  to  recover  from  the  scaling 
down  procedure  and  observe  the  pressure  at  the  shell.   The 
pressure  fluctuated  about  zero  with  an  amplitude  of 
approximately  two  atmospheres  and  analysis  of  the  positions 
of  the  groups  revealed  no  effect  of  the  initial  conditions. 
This  model  micelle  was  chosen  as  the  starting  point  for 
simulation  1. 
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Simulation  2  was  started  from  this  same  model  micelle, 
after  increasing  the  mass  of  all  head  groups  by  a  factor  of 
seven  and  allowing  10,000  time  steps  for  equilibration. 
This  mass  was  chosen  to  represent  a  common  ionic  head  group, 
the  sulfate  ion.   Only  the  mass  of  the  head  group  was 
changed;  all  of  the  intermolecular  and  intramolecular 
interactions  remained  the  same  as  in  the  previous 
simulation.   The  expected  result  of  this  would  be  simply  a 
change  in  the  dynamics  of  the  head  group.   The  program  used 
to  compute  simulations  1  and  2  is  listed  in  Appendix  D. 

In  simulation  3 ,  the  change  in  the  head  group  was  more 
extensive.   An  attempt  was  made  to  represent  the  sulfate 
head  group  in  all  of  its  intermolecular  and  intramolecular 
interactions.   The  mass  of  the  group  remained  at  seven  times 
that  of  a  methylene,  while  the  diameter  was  increased  to 
2.45  times  that  of  a  methylene.   The  equilibrium  bond  length 
and  bond  angle  for  the  head  group  were  changed,  as  were  the 
force  constants  for  bond  vibration,  bond-angle  bending,  and 
bond  rotation  involving  the  head  group.   These  values  are 
given  in  Table  6  and  compared  with  the  methylene  values. 
The  derivation  of  the  head  group  intramolecular  interaction 
parameters  based  on  the  work  of  Muller  and  coworkers  (Muller 
and  Nagarajan,  1967;  Muller  et  al.,  1968),  Cahill  et  al. 
(1968),  and  Blukis  et  al .  (1963)  is  given  in  Appendix  E. 
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Table  6 
Bond  Parameters  of  "Sulfate"  and  "Methylene"  Groups 


Parameter   "Methylene"  Value  "Sulfate"  Value 


Units 


y« 

9.25xl05 

2.7xl04 

Joule/A^/mole 

bo 

1.539 

2.6 

Angstrom 

Vb 

1.3xl05 

9.1xl05 

Joule/mole 

00 

112.15 

140 

degree 

Yr 

8313 

20000 

Joule/mole 
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Due  to  the  increase  in  the  size  of  the  head  group,  the 
range  of  head  group  interactions  was  significant  in 
comparison  to  the  size  of  the  micelle  and  truncation  of  the 
potentials  could  no  longer  be  justified.   Therefore, 
neighbor-listing  and  truncation  of  intermolecular  potentials 
were  not  employed  in  this  simulation.   It  was  also  found  to 
be  necessary  to  decrease  the  time  step  to  maintain  stability 
in  the  energetics  of  the  system.   It  is  for  these  two 
reasons  that  the  computing  efficiency  (simulation  time  per 
CPU  time)  of  simulation  3  is  dramatically  lower  than  the 
other  simulations  (Table  5) .   The  program  used  to  compute 
simulation  3  is  listed  in  Appendix  F. 

Simulation  4,  the  hydrocarbon  droplet,  was  conducted  in 
the  same  manner  as  simulation  1,  except  that  the  head  group 
was  replaced  by  a  chain  segment  (all  segments  have  the  same 
properties  to  simplify  computation)  in  its  interactions  with 
the  shell  and  other  groups.   The  radius  was  scaled  down  to 
11.928  Angstroms  to  achieve  the  liquid  hydrocarbon  density 
of  .7176  g/cc.   Since  there  was  no  sorting  out  of  the  head 
groups  for  separate  interactions,  the  efficiency  of  this 
simulation  was  18  percent  greater  than  that  of  the  first  two 
simulations. 

The  simulation  computations  were  conducted  on  a  Control 
Data  Corporation  (CDC)  Cyber  205  computer  at  the 
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Supercomputer  Research  Institute  (SCRI)  in  Tallahassee, 
Florida.   Access  to  the  SCRI  facilities  was  provided  through 
a  grant  from  the  United  States  Department  of  Energy. 

All  of  the  simulation  programs  were  modified  for  the  CDC 
FORTRAN  200  Vector-optimizing  compiler.   By  employing 
parallel  processing  techniques  when  possible,  this  compiler 
provides  very  efficient  operation  from  code  which  is  highly 
compatible  with  ANSI  standard  FORTRAN.   One  important 
difference  lies  in  the  default  word  length  of  the  Cyber  205. 
Real  numbers  on  this  computer  are  eight  bytes  in  length, 
compared  to  four  bytes  on  most  other  computers.   To  achieve 
eight-byte  precision  in  standard  FORTRAN,  the  REAL*8 
variable  type  is  used.   The  programs  in  Appendices  D  and  F 
achieve  this  level  of  precision  on  the  Cyber  205  with 
default  variable-type  coding. 

Further  detail  about  the  simulations,  along  with  the 
results  of  their  analysis,  is  given  in  the  next  chapter. 


CHAPTER  5 
RESULTS  OF  MOLECULAR  DYNAMICS  SIMULATION 

In  this  chapter  the  results  of  structural  and  dynamic 

analyses  performed  on  the  output  of  the  simulations  are 

reported.   Certain  of  these  results  are  reported  as 

time-averaged  properties;  others  are  given  as  the  change  in 

a  property  with  time.   The  former  are  averaged  over  the 

entire  length  of  the  equilibrium  run  (Table  5) ;  the  latter 

appear  in  the  figures  over  a  common  time  period  for  ease  of 

comparison.   The  simulations  are  referenced  by  number.  Run  1 

being  the  case  with  head  groups  of  the  size  and  mass  of  a 

chain  segment  (methylene) ;  Run  2  the  case  with  head  groups 

of  the  same  size  as  a  chain  segment,  but  a  mass  seven  times 

greater;  Run  3  the  case  with  head  groups  of  2.45  times  the 

size  and  seven  times  the  mass  of  a  chain  segment,  and  a 

revised  head  group  intramolecular  potential;  and  in  Run  4 

the  head  groups  replaced  by  chain  segments  to  model  a 

hydrocarbon  droplet. 

5.1  Mean  Radial  Positions  of  Groups 
The  measurement  of  mean  radial  positions  of  the  groups 
on  the  surfactant  chains  is  one  of  the  simplest,  yet  most 
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telling,  evaluations  of  micelle  structure  that  can  be  made 
on  simulation  output.   While  this  measurement  cannot  be 
obtained  reliably  for  all  groups  by  experiment,  the  position 
of  groups  is  a  fundamental  attribute  of  structural  models  of 
micelles  (Gruen,  1981) .   The  mean  radial  position  for  group 
i  (i=l  for  head,  i=9  for  tail)  is  given  by 

The  quantity  Nj^(r)  is  the  number  of  occurrences  of  the 
center  of  a  group  i  at  a  radial  position  (relative  to  the 
center  of  mass  of  the  micelle)  r  at  time  t  and  the  angle 
brackets  denote  the  average  over  the  time  period  of  the 
simulation.   The  unsubscripted  N  is  the  number  of  molecules 
in  the  system,  equal  to  the  sum  of  Nj^(r)  over  all  values  of 
r  and  i.   The  standard  deviation  corresponding  to  this  mean 
is  given  by 


o . 


-j^    <NXr)>{r-Rydr 


(5.2) 


The  mean  radial  positions  of  the  centers  of  the  nine 
groups  bracketed  by  one  standard  deviation  in  each  radial 
direction  are  plotted  for  the  four  simulations  in  Figures 
5-1  through  5-4 .   Assuming  a  normal  distribution  of  a 
group's  position  over  time,  the  range  shown  for  each  group 
in  the  figures  includes  68  percent  of  its  positions,  while 
twice  this  range  includes  9  5  percent. 
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Although  no  group  has  its  mean  position  in  the  center  of 
the  micelle,  this  does  not  indicate  a  void  within  the 
micelle.   Rather,  it  is  an  artifact  of  the  coordinate  system 
chosen.   Each  chain  segment  has  a  diameter  of  3.56 
Angstroms,  and  subsequently  excludes  a  spherical  volume  of 
this  diameter  to  any  other  group.   Since  the  spherical 
coordinate  system  provides  an  available  volume  which 
decreases  with  r,  the  excluded  volume  effect  results  in  low 
values  of  N-^Cr)  near  the  center  of  the  micelle.   Thus,  no 
group  has  its  average  position  near  the  center,  but  the 
center  is  always  within  the  excluded  volume  of  one  of 
several  different  groups.   A  range  of  three  standard 
deviations  about  the  mean  includes  approximately  all  of  a 
group's  positions,  and,  in  Figures  5-1  through  5-4, 
approaches  within  a  segment  diameter  of  the  center  for 
several  chain  segments  in  each  simulation. 

The  mean  position  results  for  Runs  1  and  2  are  virtually 
identical.   This  is  to  be  expected,  since  the  two  models 
differ  only  in  the  mass  of  the  head  group,  a  difference  that 
should  manifest  itself  in  the  dynamic  behavior  of  the 
system,  but  not  in  a  time-averaged  result.   In  these  two 
simulations,  the  head  groups  are  found  at  the  surface  of  the 
micelle  with  a  relatively  small  deviation  from  the  mean. 
Progressing  along  the  chain  toward  the  tail,  the  mean  radial 
positions  decrease  and  the  deviations  remain  roughly 
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constant  through  group  five.   The  mean  positions  of  the  last 
four  groups  are  at  approximately  the  same  radius,  but  toward 
the  tail  the  deviations  increase,  a  consequence  of  the 
excluded  volume  effect  near  the  center  of  the  micelle. 

The  results  for  Run  3  are  markedly  different.   Figure 
5-3  shows  larger  values  of  mean  positions  for  the  chain 
segments  and  greater  deviations  for  the  head  group  and  its 
adjacent  three  groups.   This  simulation  differs  from  Run  2 
in  the  size  of  the  head  group  and  the  nature  of  the 
head-chain  bond.   With  a  diameter  of  8.7  2  Angstroms,  the 
larger  head  group  creates  a  steric  effect  which  brings  about 
greater  disorder  in  the  micelle.   The  head  groups  are  spread 
over  a  greater  range  of  radial  positions,  altering  the 
positions  of  the  chain  segments  from  those  observed  in  Runs 
1  and  2 . 

Figure  5-4  shows  the  results  for  Run  4,  the  hydrocarbon 
droplet.   The  mean  positions  and  deviations  are  symmetric 
about  group  5,  the  center  of  the  molecule.   This  is  the 
proper  result,  since  the  two  ends  of  a  chain  are 
indistinguishable  from  one  another.   The  chain  ends  have  a 
slightly  larger  mean  position  and  deviation,  a  result  of  the 
greater  mobility  of  the  chain  end  relative  to  the  interior 
groups  of  the  chain.   Overall,  the  means  and  deviations  are 
quite  uniform,  indicating  a  random  arrangement  of  the 
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molecules  within  the  droplet.   Again,  the  excluded  volume 
effect  places  the  mean  positions  of  all  groups  over  two 
segment  diameters  away  from  the  center  of  the  droplet. 

In  a  simulation  of  a  micelle  of  forty  dodecyl  surfactant 
chains,  the  same  shell  model  was  used  to  surround  the 
micelle  and  the  head  groups  were  identical  to  Run  1  in  size, 
mass,  and  potentials  (Woods  et  al.,  1986).   Comparison  with 
this  simulation  reveals  the  effect  of  chain  length  on  the 
micelle  model.   In  Figure  5-5,  the  average  radial  positions 
of  Woods'  simulation  are  shown  next  to  those  of  Run  1, 
pairing  the  groups  of  Run  1  with  the  nine  groups  furthest 
from  the  head  group  of  the  longer  chain's  thirteen.   The 
dodecyl  micelle  being  larger  than  those  of  the  present  work, 
its  average  radial  positions  are  greater.   Qualitatively, 
however,  the  trend  observed  in  the  radial  positions  of  the 
groups  and  their  standard  deviations  is  quite  similar  in  the 
two  simulations. 

In  Figure  5-6,  the  average  radial  positions  of  Run  1  are 
compared  to  those  resulting  from  a  statistical  model  of  3  0 
twelve-member  chains  with  one  end  of  each  fixed  at  the 
surface  (Gruen,  1981).   In  this  model,  all  possible 
configurations  of  the  surfactant  chains  were  sampled. 
Though  Gruen 's  micelle  is  again  larger  than  that  of  Run  1, 
the  results  of  the  two  models  show  notably  similar  radial 
position  behavior  on  a  qualitative  basis. 
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A  model  micelle  of  15  sodium  octanoate  molecules,  along 
with  surrounding  molecular  water,  was  the  subject  of  another 
molecular  dynamics  simulation  (Jonsson  et  al.,  1986).   The 
sodium  carboxylate  head  group  was  explicitly  modeled  and 
surrounding  water  molecules  were  included  in  the  simulation. 
Jonsson  found  it  necessary  to  reduce  the  charge  on  the  head 
groups  to  produce  acceptable  results.   Run  3,  with  the  model 
sulfate  head  group,  is  compared  to  Jonsson 's  reduced  charge 
simulation  in  Figure  5-7.   Accounting  for  the  difference  in 
the  size  of  the  two  micelles,  the  radial  positions  of  the 
chain  segments  compare  favorably. 

5.2  Probability  Distributions  of  Group  Positions 

The  positions  of  groups  within  the  micelle  can  be 
studied  in  more  detail  by  evaluating,  at  each  radial 
position,  the  probability  of  finding  a  given  group  at  any 
given  time.   The  probability  density  function  used  in  this 
analysis  is  the  time-averaged  number  of  occurrences  of  the 
center  of  a  group  i  within  a  spherical  volume  element 
centered  at  r: 

<yv,(r)> 

PSr)  =  - TT-  (5.3) 

Evaluating  this  function  over  all  values  of  r  generates  a 
probability  density  distribution  for  each  molecular  group. 
The  results  of  this  analysis  for  each  of  the  four 
simulations  are  given  in  Figures  5-8  through  5-11.   The 


91 


c 

0 

i~i 

|M 

^ 

c^ 

0)  T3 

C 

(B 

-H 

0 

r- 

r" 

Q) 

a) 

1*^ 

4-1 

jj 

iw 

'uT 

O 

-H 

m 

H 

c 

o 

a 

•1-1. 

^ 

J-) 

p 

^  o 


-H  <M 

tn  o 

c 

r— 

>■  s 

4-J    O 

•H  a 


c 


r-i   I— J     fd 
-H    C   -iJ 

"(5  ~  0) 

o  4J  ^ 
O  C  +J 
S-i    r3 

O  -H 

(0 
•   M-l   CTi 
CO    U 

I   p  a 

in  cfl  p 

o 

0)    0)    ^ 

•H  4-1  c 
!i(    O    (0 


Ti 


S2 


CM 

=    3 

~' 

0  o 

01  cn 

_  00 

^ 

0  -c 

_ 

-H    0) 

_  ^ 

of  the  n 
is   the   h 

^y^ I 

c 

_ 

0  a 

y^^Jfi 

•H    3 

_  ^ 

J->   o 
3    ^ 

.^  ^y     ^       ^^'■^/■y 

■" 

X2  O 

r  "^^-"^""^  v^ 

£ 

< 

S-I 

CM  ■^  V^.'-''''^     -^  li^/y          1 

'    ^           • 

V/'^s^^^^    '^  ^^ 

C/f 

"5  (N 

'^t>"S>>^ 

_  o 

3 

-H 
T!    C 

in 

^/^<^y ^/^    \^ 

-  OQ 

O 

nsity 

of   Ru 
up. 

(U  in  o 

y 

73   (U   ^^ 

"^r  ^ 

X'V^'^"  ^         ^\^ 

~ 

,H  cn 

* 

>i  3 

j-"""'^        ■^        J^^'"^--^                                  ^    "■^                \ 

-P     U   rH 

r^ 

'O  'V 

^^^          ^~^~^                             N        \ 

-  CD 

•--i    Q)  -H 

; 

^    y^^              '"^^           ^  \ 

.H    rH     (C 

z' 

^x^ 

^'~"-^\        \' 

-H    0   +J 

^ 

ro         0) 

j^'^^ 

\ 

^~-~ 

XJ  -P  £ 

/ 

' 

*'*^*-^-^                                ^^            \ 

-  -"^ 

O    C  -P 

/ 

"^^S 

— 

^-1   to 
&.  -P    W 

< 

""  -■ 

^     -                              ^^\          '^ 

0  -H 
(0 

V 

-  (N 

C\    U 

\ 

-^                    \ 

1   3  a 

> 

^             \ 

in  m  3 

J 

0 

0)    0)    S-I 

S-I  £   cn 

1 

1 

1        1 

1 

1 

1      1 

1      1      1      1      1      1 

-  O 

3  -P 

&i         T3 
•H   4-1    C 

in 

-i- 

ro 

CM        -- 

o 

O) 

CO      t^ 

CD       lO       -^i-       ro       CM       T-       O 

1 — 

^ — 

^ — 

^ —        ^ — 

1 — 

o 

o     o 

O        O        O        O        O        O        O 

Cm    0    (0 

o 

o 

O 

o     o 

o 

o 

q      o 

o      q     o     o     q      o      o 

d 

d 

d 

d     d 

d 

d 

d     d 

d     <D     d      d     d      d      I 

D 

•V  "rf 


93 


a 

c 

--    n. 

O    F 

^    0 

CP  i-. 

c? 

0) 

C  TD 

•r-i     fQ 

C    Q) 

r^ 

0) 

i:  oj 

4J  jr 

4-2 

(M 

0  m 

-^ 

Ul 

C  i-l 

0 

t-i     Cm 

4-i    n 

3   0 

J2    !-i 

•H   U 

< 

en    • 

oo 

-^  n 

2J 

'^ 

r~ 

TD 

>i  5 

D 

i!   K 

Ql 

•^ 

m  «-!    . 

c  0  a 

CJ          s 

t:  M  0 

0)    L| 

>i.-i  en 

+J   3 

■^   U  <-\ 

rH    (1)  -H 

-H  iH    (C 

J2    O  +J 

ro  E 

£        0) 

O  4J  J= 

^    C  -1-! 

cu  ro 

+J    02 

U  -H 

•  ro 

C   H-i   C\ 

rH     ^4 

1   a  a 

ID    W    P 

0 

Q)    0)    J-l 

^  ^  tn 

3  -p 

Cn      T! 

■H  (W    C 

Cn   O   fO 

•V  "(i 


94 


o 

CN 

• 

Z) 

■^ 

' 

r* 

5 

? 

0) 

cc 

(U 

K 

tH 

*~ 

> 

0 

U-l 

B 

H-J 

0 

U 

0 

'Ji 

ro 

CD 

0) 

Q) 

"*" 

4J 

iH 

r* 

^ 

-H 

<W 

u 

iH 

0 

G) 

, — i 

0) 

"^ 

cn 

C 

r- 

c 
o 

s 

+J 

•H 

c  y-i 

-P 

0 

0 

g 

X! 

CM 

5-1 

m 

■■ 

•H 

ro 

& 

• 

u 

u 

p 

o< 

0 

0 

"S 

^ 

u 

—       CO 

-r-! 

TD 

& 

c     — 

TT 

>i 

""    .— 

<-• 

r— 1 

'-\J 

>. 

ro 

•^-J 

4_; 

0) 

c 

Dl 

•p-i 

^ 

-H 

cn 

JJ 

e   • 

r- 

c 

S-l    CU 

ti) 

u-l 

<n  :3 

-C 

O 

■P    0 

u 

>  U2 

0)  cri 

-P 

a 

r-; 

CO 

-H 

F-^ 

+J    rH 

i—' 

0 

-H 

•r-i 

S-1 

w  ro 

_(2 

t? 

CJ  4-1 

fD 

■P 

£ 

0    Q} 

■<j- 

C 

Q) 

C  JZ 

5-1 

^ 

CJ  +J 

cu 

_n 

TD 

ro 

m 

x: 

C\    -r-i 

* 

m 

\ 

CN' 

rH 

-H 
3 

iH   Ci 

1 

&iiH  a 

ID 

c 

Q)    n 

•H 

X5    0 

0) 

-P 

ro  5-1 

O 

P 

D 

•rH 

nH    tji 

Cr>TS 

Q)  T3 

-H 

c 

x:  c 

fc. 

■H 

Eh  ro 

o     o 


95 

probability  density  function  is  indefinite  at  the  center  of 
the  micelle,  but  by  choosing  a  finite  value  of  ^r   the 
indefinite  result  is  avoided.   Some  scatter  is  introduced  in 
the  distributions  at  small  r,  however.   Division  by  the 
volume  element  eliminates  the  excluded  volume  effect  that 
produced  no  mean  radial  positions  near  the  center  of  the 
micelle,  as  described  in  Section  5.1.   Rather,  the 
probability  density  distributions  have  peaks  throughout  the 
micelle  for  the  various  groups. 

As  in  the  case  of  the  mean  positions,  the  results  for 
Runs  1  and  2  are  quite  similar  and  show  the  same  degree  of 
order.   The  head  groups  in  these  two  simulations  remain  near 
the  surface  of  the  micelle.   Moving  along  the  chain  toward 
the  tail,  the  peak  for  each  group  occurs  closer  to  the 
center  and  the  distribution  broadens.   There  is  some 
probability  of  finding  each  of  the  groups  at  the  surface, 
but  there  is  no  probability  of  finding  the  head  group  and 
the  adjacent  three  groups  within  a  segment  diameter  of  the 
center.   Two  molecular  configurations  can  lead  to  this 
result.   First,  with  the  head  group  at  the  surface,  the 
chain  is  directed  inward  from  the  surface  and  then  is  bent 
around  such  that  its  tail  is  at  the  surface.   Second,  a 
molecule  can  lie  along  the  surface.   The  output  shows  that 
both  of  these  configurations  do  occur  in  Runs  1  and  2. 
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The  results  of  two  analogous  simulations  of  micelles  of 
40  and  52  dodecyl  surfactants  (Woods  et  al.,  1986)  give 
distributions  which  are,  except  for  micelle  radius  and 
surfactant  chain  length,  remarkably  similar  to  the  results 
of  Runs  1  and  2 .   These  simulations  also  exhibited 
comparable  average  radial  position  behavior  in  the  previous 
section.   For  the  micelle  model  used  in  these  four 
simulations,  the  difference  in  surfactant  chain  length, 
other  than  increasing  the  size  of  the  micelle,  has  a 
negligible  effect  on  the  arrangement  of  chain  segments  in 
the  micelle  interior. 

Figure  5-10  better  displays  the  disorder  present  in  Run 
3  which  was  evidenced  by  the  mean  positions  in  Figure  5-3. 
The  distribution  of  head  group  probabilities  has  two  major 
peaks:  one  at  the  surface,  and  one  approximately  one  head 
group  diameter  from  the  surface.   Groups  two  through  four 
have  roughly  the  same  peak  arrangement,  but  the  distribu- 
tions become  more  random  toward  the  end  of  the  chain.   With 
the  larger  head  group  size,  the  excluded  volume  effect  at 
the  surface  forces  one  or  more  head  groups  into  the  micelle 
interior  and  disrupts  the  orderly  arrangement  of  chain 
segments  found  in  Runs  1  and  2 . 

Figure  5-11  shows  the  distributions  for  Run  4,  the 
hydrocarbon  droplet.   Since  the  symmetry  of  the  chains  was 
demonstrated  in  the  previous  section,  the  combined 


97 

probabilities  for  the  indistinguishable  pairs  of  groups  are 
plotted.   All  groups  have  a  finite  probability  of  being 
found  anywhere  between  the  center  and  the  surface  of  the 
droplet,  and,  within  statistical  uncertainties,  the 
distributions  indicate  a  random  arrangement  of  the 
molecules. 

While  the  probability  density  function  serves  well  to 
remove  the  effect  of  excluded  volume  at  the  center  of  the 
micelle  in  evaluating  group  positions,  the  r^  dependence 
causes  some  difficulties  in  interpreting  the  results.   The 
distribution  curves  do  not  enclose  equal  areas,  and  they  do 
not  correspond  to  the  mean  radial  positions  of  Section  5.1. 
In  that  section  it  was  assumed  that  the  distribution  of 
Nj^(r)  was  normal  in  assigning  a  percentage  of  the 
distribution  to  the  standard  deviation.   By  multiplying  the 
probability  density  function  by  r^ ,  a  distribution  is 
generated  which  will  enclose  equal  areas  for  each  group,  and 
will  serve  to  evaluate  the  assumption  of  symmetry  made  in 
the  previous  section.   It  is  apparent  from  the  probability 
densities  that  the  three  groups  adjacent  to  the  head  group 
follow  the  head  group  while  the  three  or  four  groups 
adjacent  to  the  tail  follow  the  tail.   By  plotting  the 
quantity  p,  r^  for  only  the  heads  and  tails  in  Figures  5-12 
through  5-15,  clarity  is  improved  without  excessive  loss  of 
information. 
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These  plots  for  Runs  1  (Figure  5-12)  and  2  (Figure  5-13) 
show  the  head  groups  to  have  a  very  nearly  normal 
distribution,  centered  at  the  mean  radial  position,  with  a 
range  of  less  than  three  standard  deviations  in  each 
direction  along  the  abscissa.   The  tails  have  broader 
distributions,  not  quite  as  symmetric,  but  with  a  range 
within  three  standard  deviations  in  each  direction.   The 
tails  have  a  finite  probability  of  radial  positions  well 
within  a  segment  diameter  of  the  center,  confirming  the 
statistical  assumptions  of  the  previous  section.      • 

The  p,r2  plot  for  Run  3  (Figure  5-14)  reveals  a 
symmetric  peak  for  the  head  group,  centered  just  above  the 
mean  position.   It  is  broader  than  those  of  Runs  1  and  2, 
with  a  smaller  peak  in  the  interior  of  the  micelle,  as 
discussed  above.   A  comparison  of  the  areas  of  the  two  peaks 
indicates  that,  on  the  average,  two  head  groups  can  be  found 
in  the  interior  of  the  micelle.   The  tail  distribution  is 
skewed  toward  the  outer  portion  of  the  micelle  since  fewer 
groups  can  be  accommodated  in  the  center  than  in  the  outer 
portion.   The  tail  distribution  does,  however,  show  finite 
probability  of  radial  positions  within  a  segment  diameter  of 
the  center.   The  result  for  the  chain  ends  of  Run  4  is 
similar  to  that  of  the  tail  in  Run  3.   There  is  probability 


103 

of  end  groups  occurring  throughout  the  hydrocarbon  droplet, 
with  the  curve  skewed  toward  the  outer  portion  of  the 
droplet. 

Simulations  of  a  micelle  of  15  sodium  octanoate 
molecules  in  water  have  been  carried  out  by  groups  led  by 
Jonsson  (Jonsson  et  al . ,  1986)  and  Watanabe  (Watanabe  et 
al.,  1988).   Jonsson  reduced  the  charge  on  the  head  groups, 
while  Watanabe  used  the  Ewald  method  to  handle  the  long 
range  electrostatic  effects.   The  chain  length  of  the 
surfactant  monomer  of  these  simulations  differs  from  that  of 
the  present  work  by  only  one  chain  segment.   By  scaling  the 
probability  density  function  with  respect  to  the  aggregate 
size,  the  result  from  these  previous  simulations  can  be 
compared  to  that  of  this  work.   Such  a  comparison  for  the 
tail  groups  is  made  in  Figure  5-16.   The  larger  micelle  of 
Run  3  exhibits  a  distribution  with  a  higher  mean  value  of 
the  radial  position.   It  is  not  symmetric,  but  skewed  toward 
the  surface  of  the  micelle.   The  tails  of  Run  3  are 
distributed  over  a  range  equal  to  that  of  the  other  two 
simulations. 

The  experimental  technique  of  small  angle  neutron 
scattering  (SANS)  can  be  used  to  obtain  information  which 
may  be  indicative  of  micelle  structure.   Scattering 
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measurements  were  made  on  aqueous  solutions  of  lithium 
dodecyl  sulfate  above  the  critical  micelle  concentration 
(Bendedouch  et  al . ,  1983a).   The  results  of  this  study 
showed  the  alkyl  core  to  have  more  order  than  the  liquid 
state,  yet  less  than  the  all-trans  model,  and  to  be  free  of 
water.   One  result  of  the  SANS  study  is  the  Fourier 
transform  of  the  scattering  amplitude  due  to  the 
distribution  of  methyl  tails.   This  can  be  written  as 

^^  =  -r"rp„.(0=-i'^dr  (5.4) 

.4(0)  N  Jo         ^  Q 

where  A(Q)  is  the  Fourier  transform  of  the  scattering 

amplitude 
P^el^)   is  the  probability  density  of  the  tail  group 

position 
Q   is  the  scattering  vector 
Plots  of  A(Q)/A(0)  for  the  SANS  result  and  Runs  1 
through  4  are  given  in  Figure  5-17.   Again  the  results  for 
Runs  1  and  2  confirm  nearly  identical  structure.   A  plot  of 
A(Q)/A(0)  for  a  uniform  distribution  of  tails  through  the 
micelle  nearly  coincides  with  the  result  for  Run  4,  the 
hydrocarbon  droplet.   As  the  size  of  the  aggregate 
increases,  the  A(Q)/A(0)  curve  is  shifted  to  lower  values  of 
Q.   The  order  present  in  the  aggregate  affects  the  shape  of 
the  curve,  as  seen  in  the  curves  for  the  highly-ordered  Run 
1  and  the  uniform  result  without  order.   The  SANS  result, 
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for  an  aggregate  of  approximately  78  dodecyl  chains,  falls 
off  much  more  rapidly  than  the  result  for  the  much  smaller 
micelles  of  Runs  1  and  2.   Its  shape,  however,  is  quite 
similar  to  the  curves  for  Runs  1  and  2 ,  suggesting  that  the 
structure  of  Runs  1  and  2  may  be  more  like  the  experimental 
result  than  that  of  Run  3. 

The  A(Q)/A(0)  curve  for  Run  3  falls  off  more  rapidly 
than  those  of  the  other  three  simulations,  due  to  the  larger 
size  of  the  Run  3  micelle.   Its  shape  is  very  much  like  the 
Run  4  result.   The   disorder  in  structure  of  the  Run  3 
micelle  may  be  greater  than  is  indicated  for  this 
experimental  result,  though  it  is  difficult  to  make  any  real 
conclusion  regarding  this  since  the  two  micelles  are  of  such 
drastically  different  size.   Such  a  size  difference  could 
affect  internal  structure,  and,  in  fact,  the  experimentally 
determined  aggregate  number  of  78  for  dodecyl  chains  yields 
a  nonspherical  micelle. 

5.3  Conformations  of  Chain  Molecules 

From  the  analysis  of  the  positions  of  the  groups  along 
the  surfactant  chain  within  the  micelle,  it  is  apparent  that 
the  majority  of  the  molecules  are  not  in  the  all-trans 
conformation  of  the  simulation's  initial  condition.   In 
order  to  achieve  the  distributions  of  positions  observed, 
chain  bending  involving  the  transition  of  bonds  from  the 
trans  to  the  gauche  conformation  must  occur.   The 
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conformations  of  the  chains,  in  terms  of  numbers  of  gauche 
and  trans  bonds,  gives  further  information  about  the 
internal  structure  of  the  micelle. 

In  this  analysis,  a  bond  is  considered  to  be  in  the 
trans  conformation  if  the  cosine  of  its  dihedral  angle  is 
less  than  -0.5  (Weber,  1978;  Woods  et  al . ,  1986): 

cos0<-O.5  (5.5) 

Thus,  a  trans  bond  in  this  analysis  has  a  rotational 
potential  between  zero  and  the  surrounding  relative  maxima 
(see  Figure  4-2). 

The  instantaneous  percentage  of  trans  bonds  versus  time 
for  each  of  the  four  simulations  is  given  in  Figures  5-18 
through  5-21.   With  six  dihedral  angles  per  molecule,  the 
transition  of  one  bond  is  represented  by  a  change  of  0.69 
percent  in  these  figures.   Changes  of  this  magnitude  occur 
almost  continuously  (within  two  femtoseconds) ,  while  changes 
on  the  order  of  five  percent  take  place  over  a  range  of  five 
to  ten  picoseconds.   The  fluctuations  in  trans  percentage 
have  a  mean  value  which  is  dependent  on  the  simulation 
model.   The  mean  value  for  Run  1  is  68±2,  for  Run  2,  70±3, 
for  Run  3,  74±2,  and  for  Run  4,  73±3.   In  their  simulation 
of  15  sodium  octanoate  molecules,  Watanabe  et  al.  (1988) 
reported  a  value  of  78±2  percent.   Jonsson  et  al.  (1986),  in 
a  simulation  of  the  same  system,  reported  a  value  of  50±3 
percent  for  their  reduced  charge  model .   In  the 
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simulations  of  dodecyl  surfactants  (Woods  et  al . ,  1986),  a 
mean  value  of  72±3  percent  for  aggregates  of  both  40  and  52 
surfactant  molecules. 

With  the  exception  of  the  group  led  by  Jonsson,  all  of 
the  above-mentioned  simulations  found  a  percentage  of  bonds 
in  the  trans  conformation  roughly  within  one  standard 
deviation  of  one  another.   It  is  interesting  to  note  that 
the  hydrocarbon  droplet,  even  at  liquid  density,  produced  a 
mean  trans  percentage  in  the  range  corresponding  to  the 
micelles,  indicating  that  no  significant  straightening 
occurs  when  the  head  group  is  constrained  near  the  surface. 
Whether  chain  straightening  occurs  due  to  the  spherical 
geometry  of  the  aggregate  is  not  known,  since  no  simulations 
of  hydrocarbon  molecules  of  this  length  in  the  bulk  liquid 
or  aqueous  solution  phases  have  been  reported. 

Figures  5-22  through  5-25  break  the  trans  fraction  down 
by  bond  number,  each  plot  showing  the  fraction  of  the 
twenty-four  bonds  of  that  number  which  is  in  the  trans 
conformation.   The  bond  numbering  scheme  sets  Bond  1  as  the 
bond  between  the  head  group  (or  segment  1)  and  chain  segment 
2 .   Bond  2  is  the  bond  between  segments  2  and  3 ,  and  its 
dihedral  angle  is  formed  by  groups  1  through  4.   The  data 
for  these  figures  was  sampled  at  one-tenth  the  frequency 
used  in  the  previous  four  figures,  causing  the  plots  to 
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Figure  5-22.  Fraction  of  trans  conformations  for  each  of 
the  rotational  bonds  of  Run  1.  Result  is  shown  for  every 
hundredth  time  step. 
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the  rotational  bonds  of  Run  2.  Result  is  shown  for  every 
hundredth  time  step. 
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Figure  5-24.  Fraction  of  trans  conformations  for  each  of 
the  rotational  bonds  of  Run  3.  Result  is  shown  for  every 
hundredth  time  step. 
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Figure  5-25.  Fraction  of  trans  conformations  for  each  of 
the  rotational  bonds  of  Run  4.  Result  is  shown  for  every 
hundredth  time  step. 
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appear  smoother.   One  significant  result  apparent  in  these 
figures  is  the  very  high  fraction  of  trans  bonds  near  the 
head  group  in  Run  3.   In  this  simulation  the  head  group  size 
and  bond  rotational  force  constant  were  greater  than  in  Runs 
1  and  2 . 

Another  indicator  of  micelle  structure  is  the  bond  order 
parameter,  S(r),  defined  by 

S(r)=<-(3cos^0-  l)>  (5.6) 

where  9   is  the  angle  formed  by  the  bond  vector  between 
adjacent  groups  on  a  chain  and  the  radial  vector  from  the 
aggregate  center  of  mass  to  the  center  of  the  bond.   The 
parameter  is  averaged  over  all  bonds  whose  centers  are  found 
at  a  distance  r  and  over  time.   An  average  over  completely 
random  bond  orientations  would  result  in  an  order  parameter 
of  zero. 

Plots  of  S(r)  are  given  for  the  four  simulations  in 
Figures  5-26  through  5-29.   In  all  cases  little  ordering  of 
bond  orientations  exists  in  the  intermediate  portion  of  the 
aggregate  between  the  center  and  the  surface.   In  Runs  1,  2, 
and  3,  the  micelle  simulations,  S(r)  approaches  unity  at  the 
surface,  indicating  that  the  bonds  near  the  surface  tend  to 
parallel  the  radial  vector.   This  also  occurs  at  the  center 
of  the  micelle,  where  a  group  whose  center  is  coincident 
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with  the  center  of  the  micelle  would  require  that  S(r)=l  for 
a  bond  attaching  it  to  another  segment.   The  same  behavior 
of  S(r)  was  found  by  Woods  et  al.  (1986). 

In  Run  4,  the  hydrocarbon  droplet  simulation,  S(r)  goes 
to  -0.5  at  the  surface,  indicating  that  most  bonds  at  the 
surface  are  perpendicular  to  the  radial  vector.   Similar 
behavior  to  the  other  simulations  is  observed  in  the 
interior  of  the  aggregate.   The  difference  in  the  bond 
orientations  at  the  surface  of  the  micelle  and  the 
hydrocarbon  droplet  lies  in  the  head  group  interaction  with 
the  shell.   In  the  micelle,  forces  act  to  maintain  the  head 
groups  at  the  surface  while  the  chain  segments  are  repelled 
from  the  shell.   This  tends  to  orient  the  bond  between  the 
head  group  and  its  adjacent  segment  perpendicularly  to  the 
shell.   In  the  droplet  all  of  the  chain  segments  experience 
the  same  repulsive  force  at  the  shell,  resulting  in  bonds 
which  lie  parallel  to  the  shell. 

Gruen  (1981)  has  obtained  the  bond  order  parameter  as  a 
function  of  the  carbon  number  from  his  statistical  model. 
Chevalier  and  Chachaty  (1985)  also  express  a  result  of  their 
NMR  study  of  sodium  mono-n-octyl  hydrogen  phosphate  micelles 
in  this  way.   The  bond  order  parameter  as  a  function  of 
group  number,  Sj^,  can  be  approximated  from  S(r)  by 
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S,  =  -^ (5.7) 

j  NXr)dr 

The  approximation  is  due  to  the  N-[  values  being  taken  at 
group  centers  while  the  S  values  are  taken  at  bond  centers. 
Since  the  difference  between  the  two  can  be  no  longer  than 
one-half  the  bond  length,  accuracy  does  not  suffer  greatly. 
Comparisons  between  these  two  studies  and  Run  3  are  made  in 
Figure  5-30.   The  NMR  result  was  given  as  Sj^  values  relative 
to  that  of  the  head  group  bond  without  a  value  for  the  head 
group  bond;  in  the  figure  the  same  head  group  value  as  Run  3 
was  used.   The  simulation  result  shows  considerable  radial 
orientation  of  the  head  group  bond,  with  groups  removed  from 
the  head  group  showing  no  preferential  orientation.   This 
corresponds  to  the  high  trans  fraction  observed  near  the 
head  group  in  Figure  5-24.   The  NMR  result  shows  a  decrease 
in  radial  orientation  going  from  head  to  tail,  but  more 
gradually  than  the  simulation.   The  statistical  model  does 
not  show  this  trend  in  bond  orientation  along  the  chain. 
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■*■        Chevalier  (scaled) 

■         Run  3 
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Figure  5-30.   Bond  order  parameter  by  chain _ 
segment.   Results  are  shown  for  the  statistical 
model  of  Gruen  (1981),  the  NMR  study  of 
Chevalier  and  Chachaty  (1985),  and  Run  3  of  the 
present  work. 
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5.4  Shape  Fluctuations 
The  shape  of  the  model  micelles  and  droplet  can  be 
investigated  through  the  moment  of  inertia  tensor,  I: 


/.V.V 

r.y 

K. 

ly^ 

/.y 

ly. 

K. 

^y 

U. 

(5.8) 


whose  diagonal  and  off-diagonal  elements  are  calculated  from 
the  mass  and  position  of  each  of  the  216  particles  in  the 
system  as  follows: 

/.x  =  I"^v(y'*-^)  .  v=l,2 216  (5.9) 

/^y  =  -J^m,.v,y,  ,  v=  1  ,2,...,216  (5.10) 

The  eigenvalues  of  this  tensor  are  the  principal  moments  of 

inertia,  l-^,    I2,  and  I3.   In  the  case  of  all  particles  being 

of  equal  mass,  such  as  in  Runs  1  and  4,  a  spherical  micelle 

or  droplet  would  result  in  I]^=l2=l3.   In  Runs  2  and  3,  where 

the  mass  of  the  head  group  is  different  than  that  of  a  chain 

segment,  a  spherical  micelle  with  the  head  groups  equally 

spaced  over  the  surface  would  yield  equality  of  the 

principal  moments. 

Figures  5-31  through  5-34  show  the  change  over  time  of 

the  ratio  of  the  maximum  principal  moment  to  the  minimum 

principal  moment.   Any  deviation  from  the  conditions  given 

above  for  equality  of  the  principal  moments  will  result  in 

this  ratio  being  greater  than  unity.   In  all  four  of  the 

simulations,  the  ratio  Imax/^min  never  attained  a  value  of 
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unity,  despite  the  aggregate  being  confined  in  a  spherical 

shell.   Rather,  it  fluctuates  about  a  mean  value  slightly 

greater  than  l.o. 

The  result  for  Run  1  shows  a  mean  value  of  Imax/Imin  of 

1.28,  with  an  approximate  range  of  values  from  1.1  to  1.4. 
The  large  fluctuations  are  fairly  regular,  with  a  period  of 
about  five  picoseconds.   Smaller  fluctuations  are  minimal, 
on  the  order  of  0.05  on  the  Imax/Imin  scale.   The  result  for 
Run  4  is  quite  similar  to  that  of  Run  1,  with  a  mean  of  1.25 
and  range  from  approximately  1.1  to  1.4.   There  is  not  a 
regular  frequency  with  which  fluctuations  occur,  however. 

Run  2  experienced  fluctuations  in  the  approximate  range 
of  1.2  to  1.7  with  a  mean  value  of  Imax/Imin  of  1-5.   Large 
fluctuations  occurred  with  a  period  of  about  four 
picoseconds.   Smaller  fluctuations  were  distinct  and 
regular,  occurring  with  a  period  of  less  than  one  picosecond 
and  an  amplitude  of  approximately  0.01  Imax/Imin-   The 
result  for  Run  3  produced  a  mean  value  of  Imax/Imin  of  1.35 
and  a  range  within  1.1  and  1.5.   Fluctuations  of  varying 
amplitude  occurred  with  a  period  ranging  from  one  to  four 
picoseconds. 

Woods  et  al.  (1986)  obtained  mean  values  of  Imax/Imin  of 
1.25  and  1.18,  respectively,  for  their  40  and  52  member 
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dodecyl  micelles.   Jonsson  et  al.  (1986)  obtained  a  value  of 
1.49  in  their  octanoate  simulation  using  the  reduced  charge 
model . 

The  micelle  simulations.  Runs  1  through  3,  exhibited 
greater  regularity  in  the  fluctuations  of  Imax/^min  than  did 
Run  4,  the  hydrocarbon  droplet.   The  head  group  interaction 
with  the  confining  shell  must  be  considered  largely 
responsible  for  this  difference.   Runs  2  and  3,  with  the 
head  groups  of  higher  mass,  had  higher  values  of  ImaxZ-'-min 
than  Runs  1  and  4,  in  which  all  groups  were  of  equal  mass. 
The  higher  mass  head  groups  contribute  to  greater  inequality 
among  the  principal  moments  of  inertia  if  they  are 
asymmetrically  arranged  about  the  center  of  mass  of  the 
micelle.   This  was  much  more  prevalent  in  Run  2  than  in  Run 
3.   The  smaller  diameter  of  the  head  groups  in  Run  2  allows 
a  higher  local  concentration  of  the  mass  on  the  micelle 
surface  than  does  the  larger  diameter  of  the  head  groups  of 
Run  3 . 

Neutron  scattering  experiments  conducted  by  Cabane  et 
al.  (1985)  on  sodium  decyl  sulfate  micelles  concluded  that 
the  average  shape  of  the  micelles  must  be  spherical,  with 
nonspherical  shapes  possible  as  fluctuations.   This  result 
is  consistent  with  the  results  of  the  simulations.   A 
theoretical  study  of  shape  fluctuations  in  spherical 
micelles  (Ljunggren  and  Eriksson,  1984)  estimated  shape 
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fluctuations  due  to  capillary  wave  effects  on  a  time  scale 
of  100  picoseconds.   The  simulations  of  this  work  were  not 
of  sufficient  length  to  comment  on  fluctuations  of  this  long 
a  period,  but,  as  demonstrated,  shape  fluctuations  of 
significant  proportion  did  occur  on  a  time  scale  of  one 
order  of  magnitude  less. 

5.5  Pair  Correlations  of  Groups 

An  experimental  measurement  that  can  be  used  to  test  the 
validity  of  structural  models  of  micelles  is  the  pair 
correlation  distribution  function.   This  is  a  measure  of  the 
probability  of  a  pair  of  groups  being  separated  by  a  given 
distance.   Cabane  et  al.  (1985)  have  measured  the  pair 
correlations  for  certain  groups  in  sodium  dodecyl  sulfate 
micelles  of  mean  size  74.   Using  deuterium  labeling  and 
small  angle  neutron  scattering  (SANS) ,  they  obtained  results 
for  tail-tail  and  yCH2-yCH2  correlations. 

Comparing  their  SANS  result  with  the  theoretical 
predictions  of  several  structural  models,  they  found  that 
the  models  consistently  predict  the  tails  to  be  closer  to 
each  other  and  the  yCH2  groups  to  be  further  from  each  other 
than  the  experimental  result.   The  models  predicted  a 
tendency  for  tails  to  be  concentrated  near  the  center  of  the 
micelle  and  yCH2  groups  to  remain  in  a  spherical  shell  near 
the  micelle  surface. 
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In  Figures  5-35  and  5-36,  the  tail-tail  and  yCH2-yCH2 
pair  correlation  results  from  Runs  2  and  3  are  plotted  along 
with  the  SANS  result  of  Cabane  et  al.   The  axes  are  scaled 
to  put  the  results  on  an  equivalent  basis  in  aggregate  size. 
Run  3,  with  its  larger  head  groups,  is  of  a  greater  radius 
(Table  5)  and  a  lower  number  density  than  Run  2.   The 
experimental  results  lie  between  those  of  Run  2  and  Run  3, 
with  Run  2  giving  slightly  lower  pair  separations  and  Run  3 
slightly  higher.   There  is  good  agreement  between  Run  3  and 
the  experimental  result  for  the  tail-tail  correlation.   The 
agreement  between  Run  3  and  the  experimental  result  for  the 
yCH2-yCH2  correlation  is  better  than  that  of  Run  2. 

The  SANS  results  show  that  tails  are  not  concentrated  in 
the  center  of  the  micelle,  but  rather  are  distributed 
throughout.   It  is  also  indicated  by  experiment  that  the 
yCH2  groups  are  not  strictly  in  the  surface  region,  but  are 
distributed  throughout  the  micelle.   The  simulation  results, 
particularly  Run  3,  exhibit  this  same  behavior. 


CHAPTER  6 
SUMMARY  AND  CONCLUSIONS 

A  molecular  thermodynamic  model  has  been  developed  to 
describe  the  formation  of  micelles  in  a  solution  of  multiple 
surfactant  components.   Each  micelle  of  distinct  size  and 
composition  is  treated  as  the  product  of  a  reaction  whose 
equilibrium  constant  is  the  result  of  the  total  free  energy 
of  micellization  for  the  micelle.   A  seven-step,  reversible 
process  is  employed  to  calculate  the  change  in  free  energy 
for  the  surfactant  molecules  changing  their  state  from  free 
monomers  in  solution  to  aggregated  molecules  in  micelles  of 
distributed  sizes  and  compositions.   Contributions  to  the 
total  free  energy  are  obtained  due  to  solvophobic 
interaction,  mixing,  surface  formation,  conformational 
change,  head  group  interactions  and  electrostatics.   With 
the  model  for  the  free  energy  as  a  function  of  temperature 
and  composition,  distributions  of  micelle  sizes  and 
compositions  can  be  generated  through  a  set  of  reaction 
equilibria  relationships. 

Where  possible,  the  free  energy  contributions  are 
related  to  comparable  processes  on  which  experimental 
measurements  have  been  made.  The  solvophobic  term  is 
obtained  from  hydrocarbon  solubility;  the  surface  term 
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is  related  to  the  interfacial  tension  between  hydrocarbon 
and  solvent  phases.   In  addition,  reasonable  assumptions  of 
ideal  mixing  of  hydrocarbon  chains  and  ideal  solution 
behavior  of  the  dilute  surfactant  solution  are  made. 
Aggregate  size  distributions  have  been  generated  from  the 
model  for  single-component  solutions  of  nonionic  surfactants 
of  different  chain  lengths  and  at  different  temperatures  and 
solution  concentrations. 

Three  parameters  are  required  by  the  model :  one  to 
estimate  the  entropy  change  due  to  conformational  changes  in 
the  surfactant  molecules  upon  micellization  and  two  to 
estimate  the  curvature  effects  on  the  surface  free  energy  of 
spherical  and  cylindrical  interfaces.   The  parameters  can  be 
fitted  to  the  mean  aggregate  size  at  the  critical  micelle 
concentration — only  two  independent  results.   The  consequent 
interrelation  among  the  parameters  gives  curvature  parame- 
ters which  are  linearly  dependent  on  the  conformational 
parameter.   Within  the  physically  plausible  range  of  values, 
the  conformational  parameter  has  a  range  which  gives 
identical  size  distributions  when  corresponding  values  of 
the  curvature  parameters  are  used. 

This  variation  of  the  parameters  reveals  a  need  for  more 
fundamental  information  to  establish  a  truly  predictive 
model  for  the  thermodynamics  of  micelle  formation.   While 
structure  and  geometry  of  the  micelle  are  fundamental  to  the 


139 

calculation  of  free  energy  contributions,  experimental 
evidence  is  insufficient  to  completely  define  the 
thermodynamics  of  micelle  formation. 

In  order  to  acquire  some  of  the  detailed  molecular-level 
information  necessary  for  the  modeling,  the  structure  of  the 
micelle  has  been  investigated  through  the  molecular  dynamics 
computer  simulation  of  model  micelles.   Simulations  of  three 
micelles  with  different  head  group  attributes  and  one 
hydrocarbon  droplet  were  conducted.   In  each  case,  the 
aggregate  consisted  of  24  nine-member,  linear  chain 
molecules.   The  solvent  was  not  explicitly  modeled. 
Instead,  the  aggregate  was  surrounded  by  a  spherical  shell. 
Interactions  of  the  molecules  with  the  shell  were  designed 
to  mimic  the  appropriate  interaction  of  a  chain  segment  or 
head  group  with  solvent. 

The  results  of  the  simulations  have  been  analyzed  for 
the  distributions  of  molecular  groups  within  the  aggregate, 
bond  conformations  and  orientations,  and  shape  fluctuations. 
Comparisons  have  been  made  among  the  runs  to  study  the 
effects  of  head  group  size  and  mass,  with  other  simulations 
to  study  the  effects  of  chain  length,  aggregate  size,  and 
simulation  technique,  and  with  experimental  results  to 
reveal  the  differences  and  similarities  between  the 
simulations  and  experiment. 
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An  increase  in  head  group  mass  by  a  factor  of  seven  had 
no  effect  on  the  structure  of  the  micelle.   The  results  for 
Runs  1  and  2  are  virtually  identical  for  the  mean  radial 
positions  of  groups,  probability  distributions  of  group 
positions,  percentage  of  bonds  in  the  trans  conformation, 
bond  order  parameter,  and  Fourier  transform  of  tail 
distribution.   There  was  an  effect  of  head  group  mass  on  the 
dynamic  quantities.   More  rapid  gauche/trans  transitions 
were  observed  in  the  first  dihedral  angle  of  Run  2  than  in 
Run  1.   Larger  and  more  rapid  fluctuations  were  observed  in 
the  principal  moments  of  inertia  of  Run  2  than  in  Run  1. 

The  model  sulfate  head  group  of  Run  3,  with  its  higher 
mass,  larger  diameter,  and  different  intramolecular 
potentials,  had  a  great  effect  on  the  structure.   With  the 
larger  head  group  size,  the  excluded  volume  effect  at  the 
surface  forced  an  average  of  two  head  groups  into  the 
interior  of  the  micelle,  disrupting  the  order  displayed  by 
the  chains  in  Runs  1  and  2.   While  this  result  is  entirely 
possible,  the  shell  around  the  micelle  did  not  allow  other 
possible  outcomes  of  the  excluded  volume  effect  at  the 
surface,  such  as  head  groups  moving  into  the  solvent. 
Though  the  distributions  are  similar  to  other  simulations, 
the  internal  structure  of  Run  3  may  not  be  accurate. 
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The  hydrocarbon  droplet  simulation  resulted  in  a  nearly 
uniform  distribution  of  groups  within  the  droplet,  and  a 
percentage  of  trans  bonds  equivalent  to  that  of  the  three 
micelle  simulations.   This  is  in  contrast  to  statements 
based  on  other  micelle  simulations  which  suggest  there  is  a 
higher  percentage  of  trans  bonds  than  that  found 
experimentally  in  bulk  liquid  alkanes  as  proof  of  order  in 
the  micelle.   The  present  results  may  mean  that  chain 
straightening  is  more  a  product  of  the  spherical  geometry 
than  head  group  interactions,  or  it  may  just  be  an  artifact 
of  simulation.   Simulations  of  aggregates  cannot  be  compared 
with  experimental  bulk  liquid  results  until  the  same 
surfactant  molecules  are  simulated  under  "bulk"  conditions. 

Comparisons  of  the  simulation  results  of  this  work  with 
other  simulations  were  generally  consistent.   In  spite  of 
different  surfactant  molecules,  aggregate  sizes  and 
simulation  techniques,  the  results  are  qualitatively  quite 
similar. 

Comparisons  with  experimental  results  also  show 
similarities,  but  there  is  not  complete  agreement.   Since 
the  experimental  quantities  are  normally  interpretations  of 
the  measurements  which  often  conflict  from  one  study  to 
another,  the  disagreements  may  not  reflect  adversely  on  the 
simulation  results. 
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The  key  to  the  success  of  the  thermodynamic  modeling  is 
better  information  about  the  entropy  change  in  transforming 
the  hydrocarbon  droplet  into  the  head-group-free  micelle. 
The  results  of  the  present  simulations  are  inadequate  to 
determine  chain  conformations  of  monomers  in  micelles, 
droplets,  and  bulk  solvents.   This  can  serve  as  the  basis 
for  future  efforts  on  both  of  these  projects.   Chain 
conformations  and  their  energy  and  entropy  effects  must  be 
better  understood.   Simulations  of  these  molecules  as  bulk 
liquids,  in  vacuum,  and  at  infinite  dilution  in  water  would 
be  quite  beneficial  in  understanding  this  phenomena. 

Future  micelle  simulations  should  incorporate  an 
improved  solvent  model,  with  less  restrictions  on  head  group 
movement  and  shape  fluctuations.   Since  changes  in  head 
group  size  had  such  dramatic  effects,  the  size  of  the 
terminal  methyl  tail  group  should  also  be  reconsidered. 


APPENDIX  A 
MICELLE  SIZE  AND  SHAPE 

As  the  aggregate  size  of  a  micelle  increases,  a  shape 
transition  must  occur  so  that  the  volume  of  the  additional 
monomers  can  be  accomodated  while  maintaining  a  uniform 
density  througout.   In  the  case  of  a  single  surfactant 
species,  a  micelle  can  remain  spherical  only  until  the 
radius  of  the  sphere  reaches  the  length  of  an  extended, 
all-trans,  surfactant  chain.   In  the  multicomponent  case, 
the  composition  of  the  micelle  also  plays  a  role  in 
determining  its  shape.   The  geometric  model  for  a  binary 
micelle  will  be  derived  to  demonstrate  this  effect. 

The  two  surfactant  components  will  be  designated  1  and 
2,  each  with  an  aggregate  number,  N,  an  extended  chain 
length,  1,  and  a  chain  volume,  v.   Component  2  will  have  the 
longer  chain  length  and  the  larger  chain  volume.   As  in  the 
single  component  case,  the  maximum  possible  radius  of  a 
spherical  micelle  is  the  length  of  the  longest  chain  present 
in  the  micelle,  so 

N  ^v^+  N-^v^  =  —R^     ,     R<l2  (A.l) 

o 
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When  the  radius  of  the  sphere  is  greater  than  the 
shorter  chain  length,  there  is  a  spherical  "core"  in  the 
center  of  the  micelle  which  must  be  devoid  of  component  1. 
In  order  to  form  a  sphere  of  such  a  radius,  a  sufficient 
amount  of  component  2  must  be  present  to  completely  fill  the 
core.   This  is  illustrated  in  Figure  A-1.   The  radius  of  the 
core,  r,  is  given  by 

r  =  R-li     ,  /?>/,  (A. 2) 

A  fraction,  f,  of  the  total  volume  of  component  2  in  the 
micelle  will  make  up  the  core. 

/A^2^2  =  ^(^-'.)  (A.3) 

Since  a  chain  of  component  2  extending  into  the  core  must 
pass  through  the  outer  portion  of  the  micelle,  there  is  an 
upper  limit  on  the  fraction  of  component  2's  volume  that  can 
make  up  the  core. 

/<1-|^  (A. 4) 

If  both  (A.l)  and  (A. 4)  are  satisfied,  a  spherical 
micelle  will  form.   Combining  (A.l)  and  (A.3)  yields  the 
expression  for  f: 


4« 


-{N^v^+N^v^) 


3 
-  / 


f  =  — 7^ (A.5) 

A/2i^2 
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Figure  A-1.   Schematic  of  a  binary  micelle  with  monomers  of 
different  chain  lengths.   Inner  dashed  circle  depicts  the 
micelle  core  devoid  of  the  shorter  chains. 
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A  "core  fraction  constraint"  on  spherical  micelle 
formation  is  obtained  by  substituting  (A. 4)  for  f  and 
solving  for  N^: 


4ff 
3 


i-„N2V2\    1- 


1      A  3 

+  /,   )   -   N.V; 


1  /     '^   2'^2 

yV|< ^ (A. 6) 


From  A-1  a  "radius  constraint"  is  written: 

N,< (A.7) 

The  possible  combinations  of  N^  and  N2  which  will  form  a 
spherical  micelle  are  bounded  by  these  two  constraints,  as 
shown  in  Figure  A-2.   There  are  four  regions  defined  by  the 
two  constraints: 

I)  Both  constraints  met  --  spherical  micelle. 

II)  Only  radius  constraint  met  —  nonspherical . 
Ill)  Only  core  fraction  constraint  met  —  nonspherical. 

IV)  Neither  constraint  met  —  nonspherical. 
Nonspherical  micelles  are  modeled  as  prolate  spherocylin- 
ders.   The  spherical  caps  and  their  cores  have  radii  of  R 
and  r,  as  in  the  spherical  case.   The  cylindrical  portion 
has  these  same  radii  and  length  L.   These  dimensions  are 
shown  in  Figure  A-3 . 
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Figure  A-2 .   Possible  combinations  of  shorter  chains  (1) 
longer  chains  (2)  in  a  binary  micelle.   Four  regions  are 
defined  by  the  two  composition  constraints  on  micelle 
geometry. 
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A-3 .   Side  view  cf  the  prolate  spherocylinder  model 
nary  micelle.   The  micelle  is  of  radius  R  and  length 


core  is  of  radius  r  and  length  L- 
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As  in  the  spherical  case,  equations  can  be  written  for 
the  total  volume 

N  ^v^  +  N-^v^  =  -^R^  +  nLR'^  (A. 8) 

o 

and  the  core  volume 

f  N  ^v^  =  '^[R-  I  ,Y  ^  n  l[R-  I  ,Y  (A. 9) 

The  two  volume  equations  can  be  combined  to  eliminate  L, 
resulting  in  a  polynomial  in  R: 


—  l,R^ l,R^ 

3  '     3  ' 


~l]-{\-f}N^v^-N,v, 


R^ 


+  2l^[N  ^v^  +  N^v^jR-  l^[N  ^Vi  +  N^v^]  =  0  (A. 10) 

For  the  spherical  micelle  of  region  I,  L  is  zero  and 
(A.l)  gives  the  value  of  R.   In  region  III,  the  R  is  set 
equal  to  I2  and  L  is  given  by  (A. 9).   In  regions  II  and  IV, 
f  is  set  equal  to  1  -  I2/I1  ^rid  (A.  10)  is  solved  for  its 
root  between  1^  and  12- 


APPENDIX  B 
HEAD  GROUP  INTERACTION  IN  A  BINARY  MICELLE 

The  potential  energy  due  to  the  interaction  of  polar 
head  groups  at  the  surface  of  the  micelle  can  be  calculated 
as  the  sum  of  the  potentials  of  all  possible  interacting 
pairs.   Considering  ionic  head  groups  to  be  point  charges 
and  nonionic  head  groups  to  be  dipoles,  the  types  of  pairs 
possible  are  dipole/dipole,  charge/dipole,  and 
charge/charge.   The  potentials  for  all  three  are  functions 
of  the  inter-group  distance,  and  for  the  first  two  are  also 
functions  of  the  orientations  of  the  dipole(s). 

Assuming  that  the  mean  distance  between  adjacent  head 
groups  is  the  inter-group  distance  for  every  adjacent  pair, 
head  goups  can  be  arranged  hexagonally  on  the  surface 
(Figure  B-1.)   In  this  arrangement,  each  head  group  has  six 
adjacent  groups,  each  at  a  distance  "r".   If  each  group  is 
taken  as  the  center  of  a  hexagon  and  the  adjacent  pairs 
formed  with  it  are  summed,  each  pair  will  be  counted  twice 
in  the  total  of  all  adjacent  pairs  formed  by  all  groups. 
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Figure  B-1.   Hexagonal  arrangement  of  head  groups  on  the 
micelle  surface,  with  equal  separation,  r,  between  adjacent 
groups. 
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^tot  =  ^II^^«      m  =  l,2,...,A/  (B.l) 


For  a  micelle  of  N  head  groups, 

N_   _6 

2 

n   =  1  .2 6 

The  distance  r  will  first  be  determined  for  head  groups 
on  a  spherical  surface,  as  in  a  spherical  micelle  or  the 
spherical  portions  of  a  spherocylindrical  micelle.   Taking  a 
great  circle  of  the  sphere  such  that  adjacent  head  groups 
lie  on  the  circle  (Figure  B-2) ,  R  is  the  radius  of  the 
sphere  and  6(N)  is  the  angle  formed  by  two  radii  meeting  the 
surface  at  the  adjacent  groups,  as  a  function  of  N.   From 
planar  trigonometry, 

6  =  2sl„-(^)  (B.2) 

Each  triad  of  mutually  adjacent  groups  defines  an 
equilateral  spherical  triangle  with  sides  of  arclength  R6. 
Each  group  is  a  vertex  of  six  triangles.   If  six  triangles 
are  assigned  to  each  group  and  all  groups  are  summed  each 
triangle  is  counted  three  times,  so  the  number  of  spherical 
triangles  is  2N. 

From  spherical  trigonometry,  the  area  of  each  triangle 
is: 

A  =  (i3a-rT)R^  (B.3) 

where  a  is  the  angle  at  each  vertex  of  the  equilateral 
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Figure  B-2 .   Vectors  and  angles  used  in  the  calculation  of 
dipole  -  dipole  interaction  potential.   The  magnitude  of  the 
separation  vector,  r2^2'  i^  ''^^e  separation  given  in  Figure 
B-i. 
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spherical  triangle.   From  the  Law  of  Cosines  for  sides  of 
spherical  triangles, 

cos6  -  cos^6 


cosa  = 


sin  6 


cos 


2-"-'(s) 


and  from  the  double-angle  relation. 


(B.4) 


(B.5) 


2sin 


2R 


-  - 

2\R 


(B.6) 


Combining  (B.3),  (B.5),  and  (B.6)  yields  an  expression  for 
the  area  of  a  triangle  in  terms  of  r  and  R: 


A 


3cos 


:(^r 


\ 

-n)R 


(B.7) 


The  area  of  a  triangle  is  also  equal  to  the  total  area 
divided  by  the  number  of  triangles: 


A 


2N 


(B.8) 


Eliminating  the  area  in  these  two  equations  and  solving  for 
r  gives  the  distance  between  groups  for  the  spherical 
surface: 


r  =  R 


4cos(^g)-2 


(B.9) 


The  surface  of  the  spherocylinder  can  be  separated  into 
those  of  a  sphere  of  radius  R  and  cylinder  of  radius  R  and 
length  L.   The  number  of  groups  on  the  surface  of  the 


155 


micelle  is  the  sum  of  the  numbers  of  groups  on  the  two 
surfaces: 

From  (B.7)  and  (B.8)  the  number  of  groups  on  the  spherical 
portion  is  given  by 

2n 


N 


sph 


3cos' 


-(^r 


(B.ll) 


-  n 


On  the  cylindrical  portion,  the  hexagonal  pattern  will 
have  the  groups  lying  on  circles  evenly  spaced  along  the 
length  of  the  cylinder.   In  this  arrangement, 

#  of  circles  =  21 

73^ 


#  of  groups  per  circle 


n 


'n-(s) 


SI 


so  that 


N 


2nL 


cyl 


/^rsin-t-^) 

Equation  (B.IO)  can  then  be  written  as  the  following 
nonlinear  equation  in  r, 

2ti  2n  L 


(B.12) 
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3cos 


- 1 


'-(if 


3r  si 


n 


"■'it.) 


-N 


(B.13) 


4- 
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which  can  be  solved  for  r  by  iterative  techniques.  The 
values  of  R  and  L  are  calculated  by  the  method  given  in 
Appendix  A. 

For  the  calculation  of  pair  potentials,  Weston  and 
Schwarz  (1972)  give: 

z  , 2,e^ 

F  =  — - — - —  (W    1  4^1 

^  charge/charge  r^  v  *--' *  ^    'y 

Ur  \2 

- z  I  eiJ.2Cos  02 

^  charge/dipole  ^  ~      2  (B.lo) 

ur  ^2 

where  ^2   i^  the  angle  between  the  directions  r:-^2    ^^^   ^2  • 
Murrell  et  al.  (1978)  give: 


Mr/^2  3(^,  •F,2)(^2-r 


12 


^dipole/dlpole  3  5  (B.16) 

u r  12  Ljr  \2 

Assuming  all  dipoles  to  be  radially  oriented,  all  of  the 
dot  products  are  two-dimensional,  so 


M  1  ^2/  \ 

£"dipoie/d.poie  =  7^^-icos6-3cos0,cos02J  (B.17) 


Dr^ 


where  the  angles  are  those  shown  in  Figure  B-2 .   Only 
adjacent  pairs  are  being  considered  in  this  model,  and  all 
pair  separations  are  equal  to  r,  so  the  subscripts  on  r  will 
not  be  used.   Since  6  is  known  (B.2),  the  angle  terms  in 
(B.15)  and  (B.17)  can  be  eliminated: 

_ -2, e^2 

^  charge/dipole  ^  r\  D  (d.  1  O  j 

zDRr 
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dipoie/dipole 


Dr' 


r 
'^'  2R 


(B.19) 


The  total  number  of  adjacent  pairs  formed  among  N  head 
groups  is  3N.   For  a  binary  micelle, 


X, 

N 

N^ 

^2 

N 

V 

Assuming  only  adjacent  pairs  contribute  to  the  total 
interaction  potential,  the  single  component  limit  is 

£,„^  =  3yV£„  ,  .v,=  l  (B.20) 

At  X2=0,  replacing  a  "1"  with  a  "2"  replaces  six  "1-1"  pairs 
by  six  "1-2"  pairs.   With  replacements  distributed  evenly 
over  the  surface,  this  continues  until  X2=  1/4.   The  next 
replacement  of  "1"  by  "2"  replaces  four  "i-i"  pairs  by  two 
"1-2"  pairs  and  two  "2-2"  pairs.   This  continues  until  X2= 
1/2.   Beyond  this,  two  "i-i"  and  two  "1-2"  pairs  are 
replaced  by  four  "2-2"  pairs  until  X2=  3/4.   Finally,  six 
"1-2"  pairs  are  replaced  by  six  "2-2"  pairs  until  X2=  1. 
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Q?f  tot 
dx  I 


=  6N[E,^-E,,) 


=  A^(2£■22  +  2£,2-4£• 


l  1 


A/(4F22-2£,2-2f,J 


=  6A/(£22-£"i2) 


0<X2<-     (B.21) 

-<.V2<-     (B.22) 
4    '  2     ^     ^ 


-<.V2<-     (B.23) 
2^4  ^ 

-<.V2^1      (B.24) 


Integrating  the  above  and  writing  on  a  per  monomer  basis: 


£:=Fi,(3-6.v,j  +  £,2(6xJ 


0<X2<-    (B.25) 


^  . .  I  2  -  4-^'  1  J  ^  ^  12I  1  +  2.V  ,  )  +  £-22|^2.v  ,  -  2 


=  E 


,,,  2-2.V,  l+£-,2(3-2x,j  +  £22l  4.V,-- 


-<X2<-    (B.26) 
4^2        ^ 


-<X2^-    (B.27) 
2    '  4 


=  £-,2(6-6x,)  +  £22(6.v,-3^ 


<X2<1    (B.28) 


By  substituting  in  the  appropriate  pair  potentials  from 
(B.14)  through  (B.16)  the  total  energy  of  interaction  of 
adjacent  pairs  of  a  micelle  of  binary  composition  is  be 
calculated  with  one  of  equations  (B.25)  through  (B.28.) 


APPENDIX  C 

PROGRAM  LISTING  FOR  SINGLE-COMPONENT 

NONIONIC  MICELLE  CALCULATION 


REAL* 8  XSUM,XN,ARG,X 

CHARACTER  CONT, P, DISP, WHAT, SHOW, SAV, NAME*8 
REAL  N , NAVG , SAMP (3,5), NSAM ( 5 ) 
COMMON  N,DGM/I/G,GC,SCONF/J/C1,CO 
COMMON  /K/U,T,ZZ/L/GVAL(6) 
WRITE(*,*) 'NUMBER  OF  CARBONS' 
READ ( * , * ) CN 

OPEN ( 9 , FILE= ' INVALS ' , FORM= ' FORMATTED ' , ACCESS= 
& ' SEQUENTIAL ' , STATUS= ' OLD ' ) 
READ (9,*)  T,U,CO,C1,G,GC,SCONF 
READ (9,*) IMIN,IMAX, (NSAM(J) ,J=1,5) 
5  WRITE (*, 100)  C1,CO,U,T,G,GC,SCONF,IMIN,IMAX 
READ(*, • (A) ' )WHAT 
IF(WHAT.EQ. 'D' )  THEN 

WRITE(*,*) 'Enter  CI,  CO,  U,   T: • 

READ(*,*)  C1,C0,U,T 
ELSEIF(WHAT.EQ. 'P' )  THEN 

WRITE(*,*) 'Enter  G,  GC,  SCONF: • 

READ ( * , * ) G , GC , SCONF 
ELSEIF(WHAT.EQ. 'R' )  THEN 

WRITE(*,*) 'Enter  min  and  max  values  of  N: • 

READ(*,*)  IMIN,IMAX 
ELSEIF(WHAT.EQ. 'G')  THEN 

GO  TO  6 
ENDIF 
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GO  TO  5 
6  WRITE (*,*) 'PRINT  GVALUES  TO  FILE?  (Y/N) ' 
READ(*, ' (A) ' )  P 
IF(P.EQ.'Y')  THEN 

WRITE (*,*)   'ENTER  NAME  OF  FILE  (1-8  CHAR.):' 

READ(*, ' (A) ')  NAME 

OPEN (UNIT=8 , FILE=NAME , FORM= ' FORMATTED • , ACCESS= 
& ' SEQUENTIAL ' , STATUS= ' NEW ' ) 

WRITE (8,81)'  • , ' N ' , • G2 ' , ' SC ' , • G4 ' , ' G6 ' , ' GM ' , 
&'LOG(X/N) ' , 'X/N' 
ENDIF 

WRITE(*,*) 'CALCULATING  DISTRIBUTION.  .  .' 
C0LG=ALOG(C0) 
C1LG=AL0G(C1) 
CALL  HCSOL(XEQ,CN) 
DG2=ALOG(XEQ*C0) 
X1=C1/C0 
XSUM=0 . 
L0=1 

SUM1=0.0 
SUM2=0.0 
IF(P.EQ.'Y')  THEN 

KOUNT=0 

K0UNT2=1 

KUP=0 
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DO  8  I=IMIN,IMAX 
N=I 
IF(KOUNT2.EQ.8  0)  THEN 

WRITE (8, 81) '  • , 'N' , 'G2' , 'SC , '04' , '06' , 'GM' , 
&'LOG(X/N) • , 'X/N' 
K0UNT2=1 
ENDIF 

CALL  DELTAG(DG2,CN) 
ARG=N* (-DGM+C1LG) -COLG 
IF(DABS(ARG) .GT.174.673)  THEN 
IF(ARG.GT.O.O)  THEN 
WRITE(*,85)  N,DGM 
GO  TO  11 
ELSE 

XN=0.0 
ENDIF 
ELSE 

XN=DEXP(ARG) 
ENDIF 

IF(L0.GT.5)  GO  TO  7 
IF(N.EQ.NSAN(LO) )  THEN 
SAMP(1,L0)=N 
SAMP ( 2 , LO ) =GVAL ( 6 ) 
SAMP ( 3 , LO ) =XN 
L0=L0+1 
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ENDIF 

X=XN*N 

XSUM=XSUM+X 

K0UNT=K0UNT+1 

IF(N.EQ.100.)  THEN 

K0UNT=11 

KUP=1 
ENDIF 
IF(N.EQ.1000.)  THEN 

KOUNT=101 

KUP=2 
ENDIF 
IF(KUP.EQ.O)THEN 

WRITE (8, 82) '  ' , (GVAL(M) ,M=1,6) , ARG/2 . 3  03 , XN 

KOUNT2=KOUNT2+l 
ELSEIF ( KUP . EQ . 1 . AND . KOUNT . EQ . 11 )  THEN 

WRITE (8, 82) '  ' , (GVAL(M) ,M=1,6) , ARG/2 . 3  03 , XN 

KOUNT2=KOUNT2+l 

KOUNT^l 
ELSEIF ( KUP . EQ . 2 . AND . KOUNT . EQ . 1 0 1 )  THEN 

WRITE (8, 82) •  ' , (GVAL(M) ,M=1,6) , ARG/2 . 3  03 , XN 

KOUNT2=KOUNT2+l 

K0UNT=1 
ENDIF 
SUM1=SUM1+XN 
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SUM2=SUM2+X 

IF(XN.LT.1.0D-24.AND.I.GT.10000)  GOTO  10009 
8    CONTINUE 
10009    NAVG=SUM2/SU]y[l 

CMON=XSUM*C0  +  CI 
WRITE(6,84)  C1,C0,U,T 
WRITE (6, 101)  G,GC,S  CONF , NAVG , CMON 
WRITE (*,*) 'DISPLAY  SAMPLE  GVALUES?  (Y/N) ' 
READ(*, ' (A) • )DISP 

IF(DISP.EQ. 'Y')WRITE(6,83) ( (SAMP (J, K) ,J=1,3) ,K=1,5) 
11    WRITE(*,*) 'CONTINUE?  (Y/N)' 
READ(*, ' (A) ' )  CONT 
WRITE(8,84)  C1,C0,U,T 
WRITE (8,101)  G , GC , SCONF , NAVG , CMON 
C         ENDFILE(UNIT=8) 
CLOSE (UNIT=8) 
ELSE 
DO  10  I=IMIN,IMAX 
N=I 

CALL  DELTAG(DG2,CN) 
ARG=N* (-DGM+C1LG) -COLG 
IF(DABS(ARG) .GT.174.673)  THEN 
IF(ARG.GT.O.O)  THEN 
WRITE (*, 85)  N,DGM 
GO  TO  12 
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ELSE 

XN=0.0 
ENDIF 
ELSE 

XN=DEXP(ARG) 
ENDIF 

IF(L0.GT.5)  GO  TO  9 
IF(N.EQ.NSAM(LO) )  THEN 
SAMP(1,L0)=N 
SAMP ( 2 , LO ) =GVAL ( 6 ) 
SAMP ( 3 , LO ) =XN 
L0=L0+1 
ENDIF 
9      X=XN*N 

XSUM=XSUM+X 
SUM1=SUM1+XN 
SUM2=SUM2+X 

IF(XN.LT.1.0D-24.AND.I.GT.1000)  GOTO  10000 
10    CONTINUE 
10000    NAVG=SUM2/SUM1 

CMON=XSUM*C0  +  CI 

WRITE (6, 84)  C1,C0,U,T 

WRITE (6,101)  G , GC , SCONF , NAVG , CMON 

WRITE (*,*) 'DISPLAY  SAMPLE  GVALUES?  (Y/N) • 

READ(*, • (A) •)  DISP 
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IF(DISP.EQ. 'Y' ) WRITE (6, 83) ( (SAMP(J,K) ,J=1,3) ,K=1,5) 
12    WRITE (*,*) 'CONTINUE?  (Y/N) ' 

READ(*, ' (A) ' )  CONT 
ENDIF 

IF(CONT.EQ. 'Y' )   GO  TO  5 

WRITE(*,*) 'SAVE  DATA  AND  PARAMETER  VALUES?' 
READ(*, ' (A) ' )SAV 
IF(SAV.EQ. 'Y' )  THEN 

REWIND (9) 

WRITE(9,*)T,U,CO,C1,G,GC,SCONF 

WRITE(9,*)IMIN,IMAX, (NSAM(J) ,J=1,5) 
ENDIF 
CLOSE (9) 
STOP 

81  FORMAT(A,2X,A,6X,A,6X,A,6X,A,6X,A,8X,A,5X,A,4X,A, 
&5X,A,E10.4,A,F4.1,A,F6.2) 

82  FORMAT(A,F6.0,1X,F6.2,2X,F6.3,2X,F6.3,2X,F6.5,3X, 
&F8.4,2X,F6.1,2X,E9.2) 

83  FORMAT (5X, ' N= ' , F6 . 0 , 5X, ' DGM= ' ,F7,3,5X, ' XN/N= ' ,E9.2) 

84  FORMAT (IX, ' Cl= ' , E12 . 5 , 4X, 'C0=' , F5 . 1 , 4X, 'U= ' , F4 . 1 , 4X, 
&'T=' ,F7.2) 

85  FORMAT (IX, 'FOR  N= ' , F6 . 0 , ' DG) M= ' , ElO . 3/lX, 
&' OVERFLOW  ON   ' , ' XN/XN  CALCULATION') 

86  FORMAT (IX, 'FOR  N= ' , F6 . 0 , ' DG) M= ' , ElO . 3/lX, 
&'X/N  SET  TO  0  ','TO  AVOID  UNDERFLOW') 
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100  FORMAT (IX, 'Current  values  of  Data  and  Parameters  are: ' 
&/6X, 'DATA:  Cl=' ,E12.5,2X, 'C0=' , F5 . 1 , 2X, 'U= ' , F4 . 1 , 2X, 
&'T=' ,F7.2/6X, 'PARAMETERS:  G=' ,F7.4,3X, ' GC= ' ,F7.5,3X, 
&'SCONF=' ,F8.4/6X, 'RANGE  of  N  is  from  ',13,'  to  ',16// 
&  IX, 'Change  values  of  D)ata,  P)arameters,  R)ange, 

&  or  G)o?'/) 

101  FORMAT (IX, 'G=' ,F7.4,3X, •GC=' ,F7.5/1X, 'SCONF=' , F8 . 4 
&  / IX, 'MEAN  AGGREGATE  SIZE  IS',F9.2,5X, 

&  'SURFACTANT  CONCENTRATION  IS',E13.6) 
END 
******************************************************** 

SUBROUTINE  DELTAG ( DG2 , CN) 

REAL  NTRANS 

COMMON  EN , DGM/H/R, XLM/I/G , GC , SCONF/ J/Cl , CO/K/U , T , DG6 

COMMON  /L/GVAL(6) 

CALL  SHAPE1(XLC(CN) ,VC(CN) ) 

NTRANS=4.1888*XLC(CN) **3/VC(CN) 

IF ( EN. LE. NTRANS)  THEN 

SURF=12 . 5664*R*R* (NTRANS/EN) ** (1./3 • ) 
DG4=SURF*GAM(CN)*(1.-G/R)/(1.3  81*EN*T) 
ELSE 
DG4=6.2832*R*GAM(CN)*(XLM*(1-GC/R)+2*R*(1-G/R) )/ 
&    (1.381*EN*T) 

ENDIF 

CALL  DHHEDS(D) 
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SC=-SCONF* (XLC(CN)/R) **2 
DGM=DG2+DG4+SC+DG6 
GVAL(1)=EN 
GVAL(2)=DG2 
GVAL(3)=SC 
GVAL ( 4 ) =DG4 
GVAL(5)=DG6 
GVAL(6)=DGM 
RETURN 
END 
C    ******************************************************* 

FUNCTION  XLC(CN) 
XLC=1.265*CN  +1.5 
RETURN 
END 
C    ******************************************************* 

FUNCTION  VC(CN) 
VC=26.9*CN  +27.4 
RETURN 
END 
C    ******************************************************* 

FUNCTION  GAM(CN) 
COMMON  /K/Z,T,ZZ 

GAM=1.381*(57.868*CN  +  117.99  -  (.059*CN  +  .1768)*T)/ 
&   (CN  +  2.4) 
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RETURN 
END 
******************************************************* 

SUBROUTINE  HCSOL(XEQ, CN) 

REAL  K 

REAL*8  X1(10),EXPK 

DIMENSION  A (2) ,B(2) 

COMMON  /K/Z0,T,Z1 

DATA  A/997 . 098 , 258 . 015/B/1603 . 72 , 9686 . 96/R/l . 987/ 

IF(CN.LT.10. )  THEN 

1=1 
ELSE 

1=2 
ENDIF 

K=A(I)*CN  +  B(I) 
EXPK=EXP(-K/(R*T) ) 
XI ( 1 ) =EXPK/ ( 1 . +EXPK) 
DO  60  J=2,10 

X1(J)=DEXP(-(1.-2.*X1(J-1) )*K/(R*T)+DL0G(1.-X1(J-1) ) ) 
TEST=DABS(X1(J)-X1(J-1) )/Xl(J) 
IF(TEST.LT.1.0E-8)  GO  TO  61 

60  CONTINUE 

61  XEQ=X1(J) 
RETURN 
END 
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******************************************************* 
SUBROUTINE  SHAPE 1 (XLC,VC) 
COMMON  AG,Z/H/R,XL 
DATA  B/4. 18879/ 
XL=0 . 0 

R=(AG*VC/B)**(l./3.) 
IF(R.GT.XLC)  THEN 

R=XLC 

XL=(AG*VC-B*XLC**3)/(3.14159*XLC**2) 
ENDIF 
RETURN 
END 
******************************************************* 

SUBROUTINE  HEDSEP(D) 

CHARACTER* 1  FLAG 

COMMON  AG,Z/H/R,XL 

EXTERNAL  F,DF 

IF(XL.GT.O.)  THEN 

CALL  NEWTON(F,DF,D, . 01 , DO , 50 , FLAG) 

IF(FLAG.EQ. 'N' )  GO  TO  49 
D=DO 
RETURN 
49  WRITE(6,400)  AG 
RETURN 
ENDIF 
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T=COS(1.04720*(l.+2./AG) ) 
D=SQRT( (4.*T-2. )/(T-l.) ) *R 
RETURN 
400  FORMAT (IX, 'NO  HEAD  SEP.  FOUND,  N=',F7.0) 
END 

C     ****************************************ie************** 
FUNCTION  F(D) 
COMMON  AG,XZ/H/R,XL 
DATA  E/6. 28319/ 
Z=(D/R)**2 

F=E/(3*AC0S( (2.-Z)/(4.-Z) ) -3 . 14159) +E*XL/ ( 1 . 732*D 
&  *ASIN(D/(2.*R) ) )  -  AG 
RETURN 
END 

Q  ******************************************************* 

FUNCTION  DF(D) 

COMMON/ H/R, XL 

Z=D/(2.*R) 

Z2=(D/R)**2 

DF=-3.62  7  6*XL*(1./(D*D*ASIN(Z) )  +1 ./ (2 • *R*D* (ASIN (Z) ) 
&**2*SQRT(1.-Z*Z) ) )-75.398*D/(R*R*(4.-Z2)**2*(3.*ACOS( 
&  (2.-Z2)/(4.-Z2) )-3.14159)**2*SQRT(l.-( (2 . -Z2) / (4 . -Z2 
&  ))**2)) 

RETURN 

END 
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SUBROUTINE  DHHEDS(D) 
COMMON  /H/R,XL/K/U,T,DG6 
DI (T)=(.0007469*T-. 8063294 )*T+2 52. 422 
CALL  HEDSEP(D) 

DG6=3*U*U* (1 .+ (D/ (2 . *R) ) **2)/ (DI (T) *D**3*T*1 . 381E-2  3 ) 
&   *9.9907E-20 
RETURN 
END 

SUBROUTINE  NEWTON ( F , DFDX , Al , EPS , XO , J , FLAG) 

CHARACTER* 1  FLAG 

FLAG= ' Y • 

X=A1 

DO  11  1=1, J 

FX=F(X) 

IF(ABS(FX) .LT.EPS)  THEN 

xo=x 

RETURN 
ENDIF 
X=X-FX/DFDX(X) 

11  CONTINUE 
FLAG= • N • 
RETURN 
END 


C  MOLECULAR  DYNAMICS  SOURCE  PROGRAM  FOR  A  SIMPLE 

C  MODEL  OF  A  MICELLE  ENCLOSED  IN  A  SHELL.  A 

C  FIFTH-ORDER  PREDICTOR-CORRECTOR  ALGORITHM  IS 

C  USED  TO  SOLVE  THE  EQUATIONS  OF  MOTION.  A 

C  SKELETAL  MODEL  DUE  TO  WEBER  IS  USED  FOR  EACH 

C  ALKANE  MOLECULE. 

C 

c 

C  INTENDED  FOR  LONGER-RANGE  POTENTIALS.   NO 

C  TRUNCATION  OF  POTENTIALS  OR  NEIGHBOR-LISTING 

C  IS  EMPLOYED. 

C 

C  ALL  HEAD  GROUP  ATTRIBUTES  CAN  BE  SET. 


C 


PARAMETER (MHEAD=7 , HDSGMA=2 . 45 , BLHEAD= . 65 , GLHEAD= 
&  1 . 03 103E3 , BAHEAD=- . 82904 , GAHEAD=2 . 17 184E3 , GTHEAD= 
&  47.733) 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 

C0MM0N/DER/X2 (216) , Y2 (216) ,Z2 (216) , X3 (216) ,Y3 (216) , 
&   Z3(216) ,X4(216) ,Y4(216) ,Z4(216) ,X5(216) ,Y5(216) , 
&   Z5(216) 

C0MM0N/F0R/FX(216) ,FY(216) ,FZ(216) 

COMMON/DISP/DAX(216) ,DAY(216) ,DAZ (216) ,X0L(216) , 
&   Y0L(216) ,Z0L(216) 
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&   FORCON 

COMMON/ENERGY/TOTE , TOTLJ , ETOR , EBON , EBEN , EINTR , TRFRAC 
&   , EXTPRS , FTOTWL 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB , NABTOT 
COMMON/SCAL/RSPHER,RWALL,DELRS,ISCALE 
COMMON/MASS/MASS (216) 
C 

C    FLAG,  IF  IFLAG.NE.l  USES  INTPOS  AND  RANDOM  VELOCITIES 
C      IF  IFLAG.EQ.l  XO,YO,ZO,AND  ALL  DERIVATIVES  ARE  READ 
C      FROM  PREVIOUS  RUN. 
C 

IFLAG=1 
C 

C  ===  SET  NUMBER  OF  PARTICLES  IN  PRIMARY  CELL 
NAM=9 
NM=24 
NP=NM*NAM 
PART=NP 
PART=NP 
NP1=NP-1 
NP2=NP-2 

NP22=.5*PART+.01 
C 

C  ===  SET  RELATIVE  MASS  OF  PARTICLES 
DO  20  1=1, NP, NAM 
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MASS(I)=MHEAD 
2  0     CONTINUE 
C 
C  ===  SET  RADIUS  OF  SPHERE  ON  WHICH  HEAD  GROUPS  ARE  ATTACHED 

RSPHER=  6.000 

DELRS=0. 00025 

ISCALE=0 

RWALL  =RSPHER+1.0 

REQUIL=1. 

FORCON=3  0. 
C 
C  ===  SET  VALUES  OF  PHYSICAL  CONSTANTS 

XNAM=FLOAT(NAM) 

WTMOL=XNAM* 14 . +2 . 

WT PART=WTMO L/ XNAM 

RSTAR=0 . 4 

EPS=419. 

GABB=9.2  5E7 

GAB=GABB/EPS*RSTAR*RSTAR 

THA=112 . 15*3 . 14159/18  0 . 

CT0=COS(THA) 

GTHAA=1.3E5 

GTHA=GTHAA/EPS 
C   ROTATIONAL  PARAMETER  DIVIDED  BY  10  FOR  SCALE-DOWN 

GTHEE=  8.314E3 
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IF(ISCALE.EQ.1)GTHEE=GTHEE/10. 

GTHE=GTHEE/EPS 

BB=0.1539 

BO=BB/RSTAR 

B2I=1. /BO/BO 

AVO=6.02  2  5E+2  3 

BOLZ=1.38054D-23 

EPSI=EPS/BOLZ/AVO 

THIRD=l./3. 

PI=3. 1415926535 
C 
C  ===  SET  DESIRED  FLUID  STATE  CONDITION 

TR=5.913 

T=TR*EPSI 

VOL=1.333333*PI*RSPHER**3 

DR=FLOAT (NM) /VOL 
C 
C  ===  SET  RUN  FLAGS  AND  PARAMETERS 

IFLG=-1 

KB=0 

KSAVE=10 

KSORT=10 

KWRITE=10 

MAXKB=90000 

XDIST=0.1 
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C 

C  ===  SET  TIME-STEP  AND  ITS  MULTIPLES 

DELTA=0. 000850 

DELSQ=DELTA* DELTA 

DELTSQ=.5*DELSQ 

TSTEP=SQRT (WTPART*RSTAR**2/EPS/1 . E21) *1 . D+12 

TT1=TSTEP*DELTA*1.D-12 
C 
C  ===  SET  PARAMETERS  IN  PREDICTOR-CORRECTOR  METHOD 

F02=3./16. 

F12=251./360. 

F32=ll./18. 

F42=l./6. 

F52=l./60. 
C 
C  ==-  SET  DISTANCES  FOR  POTENTIAL  CUT-OFF,  VERLET  LIST,  ETC, 

CUBE=2 . *RSPHER 

CUBE2=.5*CUBE 

RC=2 . 5 

RCHH=2 . 5 

RLIST=(RC+.25)**2 

RDMAX=CUBE2  *CUBE2 

IF(RLIST.GT.RDMAX)  RLIST=RDMAX 

IF(RC.GT.CUBE2)  RC=CUBE2 
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C  ===  SCALE  FACTOR  FOR  VELOCITIES  DURING  EQUILIBRATION 

AHEAT=DELSQ*PART*3 . *TR 
C 
C  ===  WRITE  TAPE  HEADING 

OPEN(2  0,FILE='C811RUN1' , ACCESS= ' SEQUENTIAL' , 
&  FORM= ' UNFORMATTED '  ) 
WRITE (20)  NP,NAM,NM,DR,TR,EPS,RSTAR 
WRITE (20)  CUBE , VOL , RC , RLIST , RDMAX , DELTA , TSTEP 
C 
C  ===  SHIFTED-FORCE  CONSTANTS 

RRMAX=1./RC 

RRMAX 6 =RRMAX * * 6 

RRMAX 9 =RRMAX * * 9 

ESHFT=2 1 . *RRMAX6-20 . *RRMAX9 

ESHFTA=18 . *RRMAX* (RRMAX9-RRMAX6) 

FSHFT=ESHFTA 
C    SHIFTED  FORCE  POTENTIALS  FOR  HEAD-HEAD  REPULSIVE 
C     INTERACTIONS 
C 

RCI=1./RCHH 

RCI2=RCI*RCI 

RCI3=RCI2*RCI 

RCI4=RCI3*RCI 

RCI9=RCI3*RCI3*RCI3 

RCI13=RCI4*RCI*RCI4*RCI4 
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HSHFT=RCI3* (13 . *RCI9+4 . ) 

HSHFTA=3.*RCI4*(4.*RCI9  +  1.  ) 

HFSHFT=HSHFTA 
C 
C  ===  CORRECTIONS  FOR  LONG-RANGE  INTERACTIONS 

RC3=RC**3 

RC6=RC3*RC3 

C0RE=4 . *PI*DR* ( 1 . /6 . /RC6-0 . 5/RC3 ) 

DE=CORE 
C 
C  ===  INITIALIZE  SUM  ACCUMMULATORS 

XSUM=0. 

SUME=0. 
C 
C  ===  PRINT  PARAMETERS 

WRITE(*,900) 
900    FORMAT (IHl///) 

WRITE(*,902) 
902    FORMAT (7X, 49( •*')  ) 

WRITE (*, 904) 
904    F0RMAT(7X,  '*'  ,T56,  '*' ) 

WRITE (*, 906)  NM 
906    FORMAT ( 7 X, •*• ,2X, 'MOLECULAR  DYNAMICS  FOR ',13, 
&  '  N-ALKANE  MOLECULES' ,T5 6,  '*' ) 

WRITE (*, 904) 
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WRITE(*,917)  NAM 
917    FORMAT (7X, •*' ,5X, 'WITH' ,13, '  PARTICLES  PER 
&MOLECULE' ,T56, '*' ) 
WRITE(*,904) 
WRITE(*,902) 
WRITE(*,904) 
WRITE (*, 908)  EPSI,RSTAR 
908    FORMAT ( 7X, •*' ,2X, 'EPSI/K  =  ' , F7 . 3 , T36 , ' RSTAR  =  ', 
&F7.3,T56, '*' ) 

WRITE(*,910)  TR,T 
910    FORMAT (7X, '*' ,2X, 'TR      =  ',F7.3,T36,'     T  =  ', 
&F7.3,T56, '*') 

WRITE(*,912)  DR,VOL 


912    FORMAT (7X, '*' ,2X, 'DR      =  ' , F7 . 3 , T36 , ' VOL 


—  I 


I 


&F8.3,T56, '*') 

WRITE (*, 914)  CUBE 

914    FORMAT (7X, '*' ,2X, 'CUBE  =  ' , F7 . 3 , T55 , ' * ' ) 

WRITE(*,916)  RC 

916    FORMAT (7X, '*' ,2X, 'RC  =  ' , F7 . 3 , T55 , ' * ' ) 

WRITE (*, 918)  RLIST,RDMAX 

918    FORMAT ( 7X, '*• ,2X, 'RLIST  =  ' , F7 . 3 , T36 , 'RDMAX  =  ', 
&F7.3,T56, '*') 

WRITE (*, 92  0)  DELTA 

920    FORMAT ( 7X, '*' ,2X, 'DELTA  =  ' , F9 . 5 , T56 , ' * ' ) 

WRITE(*,904) 
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WRITE(*,922)  TSTEP 
922    FORMAT  (7X,  '*•  ,2X,  'TIME  UNIT  =  ',F6.3,'E-12  SEC, 
&T56, '*•) 

WRITE(*,924)  TTl 
924    FORMAT  (7X,  '*'  ,2X,  'TIME  STEP  =  ',1PE10.3,'  SEC, 
&T56, '*• ) 

WRITE(*,904) 
WRITE (*, 92  6)  DE 
926    FORMAT ( 7X, •*• ,2X, 'ENERGY    CORRECTION  =  ',F7.3,T56, 

WRITE(*,904) 

WRITE (*, 902) 

WRITE (*, 932)  RSPHER 
932     FORMAT (////lOX, 'RADIUS  OF  SPHERE  =  ',F6.2) 

WRITE (*, 934)  GTHE 
934     FORMAT (//lOX, 'ROTATIONAL  POTENTIAL  PARAMETER  =  ', 
*F7.3/) 


C 
C 

c 
c 


OPEN (30, FILE= ' LASTDAT ' , ACCESS= ' SEQUENTIAL ' , FORM= 
&' UNFORMATTED') 

IF(IFLAG.EQ.l)  GOTO  199 


C 
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C  ===  LOAD  INITIAL  POSITIONS  OF  ATOMS 
CALL  INTPOS(RSPHER) 
CALL  POS  PRI  ( KB ,  irA ,  IT A>I ) 
C 

C  ===  PRINT  RUN-TABLE  HEADING 
V?RITE(*,S30) 
930    FORMAT (1H1////4X, 'KB' ,5X, 'RSPH' ,4X, ' ENRG ' ,5X, 'El' , 
&    4X, 'DIST' ,6X, 'TEMP' ,3X, ' TRFRC ' ,3X, 'NAB' ,3X, 
&'TOT  ENR' ,3X, 'FTOTWL' ,3X, 'PRESSURE'/) 
C 
C  ==  LOAD  INITIAL  VELOCITIES  OF  ATOMS 

CALL  INTVEL(AHEAT,PART) 
C 

C  ===  ASSIGN  INITIAL  ACCELERATIONS  BASED  ON  INITIAL 
C       POSITIONS 

CALL  EVAL(RWALL) 
C 

C  ===  SCALE  ACCELERATIONS  AND  STORE  STARTING  POSITIONS 
DO  530  1=1, NP 

X2 (I)=FX(I) *DELTSQ/MASS (I) 
Y2 ( I ) =FY ( I ) *DELTSQ/MASS ( I ) 
Z2 (I) =FZ (I) *DELTSQ/MASS (I) 
53  0    CONTINUE 
GOTO  188 
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C   READ  POSITIONS  AND  DERIVATIVES  INSTEAD  OF  USING 

C   FCC,  INVEL,  ACCELERATION  

C 

199      WRITE(*,930) 

READ (30)  KB,IFLG,NM,NAM,RSPHER,XSUM,SUME 
RWALL=RSPHER+1 . 0 
DO  200  K=1,NP 

READ(30)  XO(K) ,YO(K) ,ZO(K) 
READ(30)  X1(K) ,Y1(K) ,Z1(K) 
READ(30)  X2(K) ,Y2(K) ,Z2 (K) 
READ(30)  X3(K) ,Y3 (K) ,Z3 (K) 
READ(30)  X4(K) ,Y4(K) ,Z4(K) 
READ(30)  X5(K) ,Y5(K) ,Z5(K) 
READ (30)  DAX(K) ,DAY(K) ,DAZ(K) 
2  00      CONTINUE 
188      DO  377  1=1, NP 
XOL(I)=XO(I) 
YOL(I)=YO(I) 
ZOL(I)=ZO(I) 
377         CONTINUE 
C 
C 

c 

C  ===  ENTER  MAIN  LOOP  OF  SIMULATION 
IF(IFLAG.NE.l)  GOTO  777 
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CALL  PREDCT(NP) 

CALL  EVAL(RWALL) 

CALL  CORR(DELTSQ) 

IF ( ISCALE . EQ . 1 )  RSPHER=RSPHER+DELRS 

RWALL=RSPHER+1 . 0 
777       NS=KB+1 

DO  599  NTIMES=NS,iyiAXKB 

KB=KB+1 

CALL  PREDCT(NP) 

CALL  EVAL(RWALL) 

CALL  CORR(DELTSQ) 
C 
C  ===  CALCULATE  MEAN  SQUARE  DISPLACEMENT  &  KINETIC  ENERGY 

SUMVEL=0. 

TDIST=0. 

DO  540  1=1, NP 

TDIST=TDIST+DAX(I) **2+DAY(I) **2+DAZ (I) **2 
SUMVEL=SUMVEL+(X1(I) **2+Yl(I) **2+Zl(I) **2) 
&*MASS(I) 
54  0      CONTINUE 

TDIST=TDIST/PART 

EK=SUMVEL/ ( 2 . *PART*DELSQ) 
C 
C  ===  ACCUMMULATE  SUMS  FOR  PROPERTY  AVERAGES 

XSUM=XSUM+SUMVEL 
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SUME=SUME+TOTE 
C 

C  ===  PROPERTY  CALCULATION  &  PRINT-OUT  AT  INTERVALS 
IF(MOD(KB,KWRITE) .NE.O)  GOTO  550 
FKB=FLOAT (KB) *PART 
TMP=XSUM/ ( 3 . *DELSQ*FKB) 
ENR= ( SUME/ FKB+CORE ) 
E 1=T0TE/ PART+ CORE 
ET0T=E1+EK 

RLTIM=DELTA* FLOAT (KB) *TSTEP 

WRITE ( * , 940)  KB, RSPHER, ENR, El , TDIST, TMP, TRFRAC , 
&  NABTOT,ETOT,FTOTWL,EXTPRS 
C 

IF(MOD(KB, 1000) .NE.O)  GOTO  550 
CALL  POSPRI(KB,NM,NAM) 
WRITE(*,930) 
940        F0RMAT(1H  , 15 , 4F8 . 3 , FIO . 3 , IX, F6 . 3 , 16 , 2X, F8 . 4 , 
&  2 (3X,F9.3) ) 
C 

C  ===  DURING  FIRST  OF  RUN,  SCALE  VELOCITIES  FOR  TEMPERATURE 
550      IF(IFLG.LT.l)  CALL  EQBRAT (SUMVEL,AHEAT, TDIST, XDIST 
&,NP,IFLG,KB,NAM) 
C 

C  WRITE  DATA  ONTO  TAPE  FOR  LATER  USE  

IF(KB.EQ.O)  GOTO  777 
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IF(IFLG.EQ.O)  GOTO  588 
IF(MOD(KB,10) .NE.O)  GOTO  588 

WRITE (20)  KB , TMP , ETOT , El , EK , TOTE , EINTR , TOTLJ , ETOR , 
&EBON,EBEN 

WRITE(20)  XO,YO,ZO 

WRITE(20)  X1,Y1,Z1 

WRITE(20)  X2,Y2,Z2 

C         WRITE(20)  X3,Y3,Z3 

C         WRITE(20)  X4,Y4,Z4 

C         WRITE(20)  X5,Y5,Z5 

588      IF(MOD(NTIMES, 100) .NE.O)  GOTO  599 
REWIND  3  0 
WRITE (30)  KB,IFLG,NM,NAM,RSPHER,XSUM,SUME 
DO  598  K=1,NP 

WRITE(30)  XO(K) ,YO(K) ,ZO(K) 
WRITE(30)  X1(K) ,Y1(K) ,Z1(K) 
WRITE(30)  X2(K) ,Y2(K) ,Z2(K) 
WRITE(30)  X3(K) ,Y3(K) ,Z3(K) 
WRITE(30)  X4 (K) ,Y4(K) ,Z4(K) 
WRITE(30)  X5(K) ,Y5(K) ,Z5(K) 
WRITE (30)  DAX(K) ,DAY(K) ,DAZ(K) 

598  CONTINUE 

599  CONTINUE 
C        CLOSE (20) 

CLOSE (30) 
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STOP 

END 

SUBROUTINE  POSPRI (KB, NM, NAM) 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

1 

11=0 

WRITE(*,942)  KB 
942      FORMAT (1H1///7X, 'POSITIONS  OF  GROUPS  AT  TIME-STEP: 
&',I6//) 

DO  4  04  JJ=1,NM 
WRITE(*,935) 
935        F0RMAT(//2X, 'MOLECULE' ,7X, 'X' ,11X, 'Y' ,12X, 'Z' , 
&12X, 'R' , 

&  9X, 'BOND  LEN'/) 

DO  4  04  KK=1,NAM 
11=11+1 

RR=SQRT(X0(II)**2+Y0(II)**2+Z0(II)**2) 
XXX=XO (II) -XO (II-l) 
YYY=YO (II) -YO (II-l) 
ZZZ=Z0(II)-Z0(II-1) 
BB=SQRT(XXX*XXX+YYY*YYY+ZZZ*ZZZ) 
IF(KK.EQ.l)  WRITE(*,909)  JJ, XO (II) , YO (II) , 
&ZO(II) ,RR 

IF(KK.GT.l)  WRITE(*,911)  XO (II) , YO (II) , ZO (II) , 
&RR,BB 
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4  04      CONTINUE 

909      FORMAT(5X,I2,4(5X,F8.4) ) 

911      FORMAT(7X,5(5X,F8.4) ) 

RETURN 

END 


SUBROUTINE  INTPOS (RSPHER) 
C 

C   SET-UP  INITIAL  POSITIONS  OF  PARTICLES.  HEAD  GROUPS  ARE 
C     ASSIGNED  TO  FIXED  POSITIONS  ON  THE  SURFACE  OF  A  SPHERE 
C     OF  RADIUS  "RSPHER". 
C 

COMMON/ POS/XO (216) ,Y0(216) ,Z0(216) 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB , NABTOT 
COMMON/ CONT/ GAB , CTO , GTHA , GTHE , BO , B2 I , Al , AG , EPS 
DIMENSION  THETA(50) ,PHI(50) 
C 

PI=3. 1415926535 
PI180=PI/180. 
C 

C   TWENTY-FOUR  MOLECULES 
THETA(1)=0. 
PHI(1)=0. 
THETA(24)=PI 
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PHI(24)=0. 

PH24=0. 

DO  70  1=2,5 

PHI(I)=PH24 

PH24=PH24+PI/2. 

70  THETA(I)=36.*PI180 
DELPHI=2.*PI/7 
PH24=DELPHI/2. 

DO  71  1=6,12 
PHI(I)=PH24 
PH24=PH24+DELPHI 

71  THETA(I)=72.*PI180 
PH24=0. 

DO  72  1=13,19 
PHI(I)=PH24 
PH24=PH24+DELPHI 

72  THETA(I)=108.*PI180 
PH24=PI/4. 

DO  73  1=20,23 
PHI(I)=PH24 
PH24=PH24+PI/2. 

73  THETA(I)=144.*PI180 
C 

60    WRITE(*,909) 
909     FORMAT (1H1////7X, 'INITIAL  POSITIONS  OF  HEAD  GROUPS' 
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&    //7X, 'MOLECULE' ,5X, 'THETA' ,8X, 'PHI'/) 
DO  101  1=1, NM 

TTT=THETA ( I ) / PIl 8  0 
PPP=PHI(I)/PI18  0 
WRITE (*, 900)  I,TTT,PPP 
900       FORMAT(10X,I3,2(6X,F7.3) ) 
101     CONTINUE 
C 

c 

C    STANDARD  DISTANCES  ALONG  TRANS  CHAIN 
XLEN=0.5*B0*SQRT(2.*(1.-CT0) ) 
ZLEN=SQRT(BO*BO-XLEN*XLEN) 
ZLENSQ=ZLEN*ZLEN 
C 

C  ===   ADD  ATOMS  TO  FORM  MOL  IN  TRANS  POS  ==== 
DO  12  0  1=1, NM 
IADD=(I-1)*NAM 
THR=THETA ( I ) 
PHIR=PHI(I) 
SR=SIN(THR) 
CR=COS(THR) 
SPR=SIN(PHIR) 
CPR=COS(PHIR) 
C 
C   ODD  NUMBERED  GROUPS  ON  A  MOLECULE 
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DO  92  J=1,NAM,2 
K=IADD+J 

RK=RSPHER-FL0AT(J-1) *XLEN 
XO(K)=RK*SR*CPR 
YO(K)=RK*SR*SPR 
ZO(K)=RK*CR 
92       CONTINUE 
C 

C   EVEN  NUMBERED  GROUPS  ON  A  MOLECULE 
DO  95  J=2,NAM,2 
K=IADD+J 

RDIS=RSPHER-FL0AT(J-1) *XLEN 
RK=SQRT (RDIS*RDIS+ZLENSQ) 
BETA=ASIN ( ZLEN/RK) 
THE=THR-BETA 
ST=SIN(THE) 
XO(K)=RK*ST*CPR 
YO(K)=RK*ST*SPR 
ZO(K)=RK*COS(THE) 
95       CONTINUE 
12  0     CONTINUE 
RETURN 
END 
SUBROUTINE  INTVEL ( AHEAT , PART) 
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C  ===  ASSIGN  INITIAL  VELOCITIES  TO  ATOMS 
C 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB , NABTOT 
COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 
COMMON/MASS/MASS (216) 
SUMX=0 . 
SUMY=0 . 
SUMZ=0. 
DO  200  1=1, NM 
XX=RANF(DUM) 
YY=RANF(DUM) 
ZZ=RANF(DUM) 
NX=(I-1)*NAM 
DO  200  J=1,NAM 
X1(J+NX)=XX 
Y1(J+NX)=YY 
Z1(J+NX)=ZZ 

SUMX=SUMX+X1 ( J+NX) *MASS ( J+NX) 
SUMY=SUMY+Y1 (J+NX) *MASS (J+NX) 
SUMZ=SUMZ+Z1 (J+NX) *MASS (J+NX) 
2  00    CONTINUE 
C 
C  ===  SCALE  VELOCITIES  SO  THAT  TOTAL  MOMENTUM=ZERO 

x=o. 

DO  210  1=1, NP 
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X1(I)=X1(I) -SUMX/PART/MASS(I) 
Yl ( I ) =Y1 ( I ) -SUMY/PART/MASS ( I ) 
Z1(I)=Z1(I)-SUMZ/PART/MASS(I) 
X=X+(X1(I)**2+Y1(I)**2+Z1(I)**2)*MASS(I) 
210    CONTINUE 
C 

C  ===  SCALE  VELOCITIES  TO  DESIRED  TEMPERATURE 
HEAT=SQRT (AHEAT/X) 
DO  220  1=1, NP 

X1(I)=X1(I)*HEAT 
Y1(I)=Y1(I)*HEAT 
Z1(I)=Z1(I)*HEAT 
220    CONTINUE 
RETURN 
END 


SUBROUTINE  PREDCT(NP) 
C 

C  ===  USE  TAYLOR  SERIES  TO  PREDICT  POSITIONS  &  THEIR 
C        DERIVATIVES  AT  NEXT  TIME-STEP 
C 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 

COMMON/ DER/X2 (216) ,Y2(216) ,Z2(216) ,X3(216) ,Y3(216) , 
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&  Z3(216) ,X4(216) ,Y4(216) ,Z4(216) ,X5(216) ,Y5(216) , 
&Z5(216) 

COMMON/FOR/FX(216) , FY (216) , FZ (216) 
DO  300  1=1, NP 


300 


X0(I 

Y0(I 

ZO(I 

X1(I 

Y1(I 

Z1(I 

X2  (I 

Y2(I 

Z2(I 

X3  (I 

Y3(I 

Z3(I 

X4(I 

Y4(I 

Z4  (I 

FX(I 

FY  (I 

FZ(I 
CONTINUE 
RETURN 
END 


=X0 
=Y0 
=Z0 
=X1 
=Y1 
=Z1 
=X2 
=Y2 
=Z2 
=X3 
=Y3 
=Z3 
=X4 
=Y4 
=Z4 
=  0. 
=  0. 
=0. 


(I)+X1(I)+X2 (I)+X3 (I)+X4 (I)+X5(I) 

(I)+Y1(I)+Y2(I)+Y3(I)+Y4(I)+Y5(I) 

(I)+Z1(I)+Z2(I)+Z3(I)+Z4(I)+Z5(I) 

(I)+2.*X2(I)+3.*X3(I)+4.*X4(I)+5.*X5(I) 

(I)+2.*Y2 (I)+3.*Y3 (I)+4.*Y4(I)+5.*Y5(I) 

(I)+2.*Z2 (I)+3.*Z3 (I)+4.*Z4(I)+5.*Z5(I) 

(I)+3.*X3(I)+6.*X4(I)+10.*X5(I) 

(I)+3.*Y3 (I)+6.*Y4 (I)+10.*Y5(I) 

(I)+3.*Z3(I)+6.*Z4(I)+10.*Z5(I) 

(I)+4.*X4(I)+10.*X5(I) 

(I)+4.*Y4 (I)+10.*Y5(I) 

(I)+4.*Z4(I)+10.*Z5(I) 

(I)+5.*X5(I) 

(I)+5.*Y5(I) 

(I)+5.*Z5(I) 
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SUBROUTINE  CORR(DELTSQ) 
C 

C  ===  CORRECT  PREDICTED  POSITIONS  AND  THEIR  DERIVATIVES 
C 

COiyiMON/POS/X0(216)  ,Y0(216)  ,Z0(216) 
COMMON/VEL/Xl(216) , Yl (216) , Zl (216) 

COMMON/DER/X2(216) , Y2 (216) , Z2 (216) , X3 (216) , Y3 (216) , 
&  Z3(216) ,X4(216) ,Y4(216) ,Z4(216) ,X5(216) ,Y5(216) , 
&  Z5(216) 

COMMON/FOR/FX(216) , FY (216) , FZ (216) 

COMMON/DISP/DAX(216) ,DAY(216) ,DAZ(216) ,X0L(216) , 
&  Y0L(216) ,Z0L(216) 

COMMON/ PARM/FO 2 , F12 , F32 , F42 , F52 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB , NABTOT 
COMMON/ PRO P/SUME , XSUM 

COMMON/SCAL/RSPHER, RWALL, DELRS , ISCALE 
COMMON/MASS/MASS (216) 
C 

IF(MOD(KB, 10) .NE. 0. OR. ISCALE. NE.l)  GOTO  610 
C        IF(RSPHER  .LE.  2.800)  GO  TO  610 
RS  PHER=RS  PHER- DELRS 
RWALL  =RWALL  -DELRS 
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610    DO  690  1=1, NP 

XERR=X2 ( I ) -DELTSQ*FX ( I ) /MASS ( I ) 

YERR=Y2 ( I ) -DELTSQ*FY ( I ) /MASS ( I ) 

ZERR=Z2 (I)-DELTSQ*FZ(I)/MASS(I) 

XO ( I ) =X0 ( I ) -XERR*F02 

X1(I)=X1(I)-XERR*F12 

X2(I)=X2(I)-XERR 

X3 (I)=X3 (I) -XERR*F32 

X4 ( I ) =X4 ( I ) -XERR*F4  2 

X5 ( I ) =X5 ( I ) -XERR*F5  2 

Y0(I)=Y0(I)-YERR*F02 

Y1(I)=Y1(I) -YERR*F12 

Y2(I)=Y2(I)-YERR 

Y3 (I) =Y3 (I) -YERR*F32 

Y4 ( I ) =Y4 ( I ) -YERR*F4  2 

Y5 (I) =Y5 (I) -YERR*F52 

Z0(I)=Z0(I)-ZERR*F02 

Z1(I)=Z1(I)-ZERR*F12 

Z2(I)=Z2 (I)-ZERR 

Z3(I)=Z3(I)-ZERR*F32 

Z4 (I)=Z4 (I)-ZERR*F4  2 

Z5 (I)=Z5 (I) -ZERR*F52 
69  0   CONTINUE 

IF(MOD(KB,KSORT) .NE.O.OR.ISCALE.NE.l)  GOTO  680 
DO  691  1=1, NP 
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RSQ=XO(I)*XO(I)+YO(I)*YO(I)+ZO(I)*ZO(I) 
R=SQRT(RSQ) 
RFACT  =  1.0-DELRS/R 
XO(I)  =  XO(I)*RFACT 
YO(I)  =  YO(I) *RFACT 
ZO(I)  =  ZO(I)*RFACT 
691    CONTINUE 
C 

C  ===  DISPLACEMENTS 
680   DO  692  1=1, NP 

DAX ( I ) =DAX ( I ) -XO ( I ) +XOL ( I ) 
DAY ( I ) =DAY ( I ) -YO ( I ) +YOL ( I ) 
DAZ(I)=DAZ(I)-ZO(I)+ZOL(I) 
C 

C  ===  STORE  NEW  POSITIONS 
XOL(I)=XO(I) 
YOL(I)=YO(I) 
ZOL(I)=ZO(I) 
692    CONTINUE 
RETURN 
END 


SUBROUTINE  EQBRAT ( SUMVEL , AHEAT , TDIST , XDIST , NP , 
&    IFLG,KB,NAM) 
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C 

C  ===  SCALE  VELOCITIES  DURING  INITIAL  TIME-STEPS 
C 

COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 
COMMON/ PRO P/SUME , XSUM 

COMMON/SCAL/RSPHER,RWALL,DELRS,ISCALE 
C 

IF(IFLG.EQ.-l)  GOTO  720 

IF(TDIST.GT.XDIST.0R.IFLG.EQ.-2)  GOTO  750 
720    HEAT=SQRT(AHEAT/SUMVEL) 
DO  730  1=1, NP 

X1(I)=X1(I)*HEAT 
Y1(I)=Y1(I)*HEAT 
Z1(I)=Z1(I) *HEAT 
730    CONTINUE 
RETURN 
C 

C  ===  AT  END  OF  EQUILIBRATION  STAGE,  SET  PROPERTY  SUMS 
C       TO  ZERO 
750    IFLG=  1 
KB=0 
SUME=0. 
XSUM=0 . 
RETURN 
END 
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SUBROUTINE  EVAL(RWALL) 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

COMMON/FOR/FX(216) , FY (216) , FZ (216) 

COMMON/NABLST/LIST( 12000) ,NABORS(216) 

COMMON/ CONT/ GAB , CTO , GTHA , GTHE , BO , B2 I , Al , AO , EPS 

COMMON/PROP/SUME , XSUM 

COMMON/ PRO PA/ RC , RLIST , FSHFT , ESHFT , ESHFTA , CUBE , CUBE2 , 
&  RDMAX , RCHH , HSHFT , HSHFTA , HFSHFT , REQUIL , FORCON 

COMMON/ENERGY/TOTE , TOTLJ , ETOR , EBON , EBEN , EINTR, TRFRAC , 
&  EXTPRS , FTOTWL 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB , NABTOT 

DIMENSION  B(216) ,BX(216) ,BY(216) ,BZ(216) 


C 

C 

INITIALIZE  ACCUMULATORS 

NABTOT=0 

ICOUNT=0 

K=0 

ITRAN  =0 

TOTE=0.0 

TOTLJ=0.0 

FTOTWL=0 . 0 

ETOR=0 . 0 
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EBON=0.0 
EBEN=0.0 
EINTR=0.0 
DO  150  1=1, NP 
FX(I)=0.0 
FY (I) =0.0 
FZ(I)=0.0 
150   CONTINUE 
C 

C   CALCULATE  LJ  FORCES 
LCHK=2 

IF(MOD(KB,KSORT) .EQ.O)  LCHK=1 
C 

C   OUTER  LOOP  OVER  PARTICLES 
DO  300  1=1, NPl 
ICHK=I-l+NA]yi 
XI=XO(I) 
YI=YO(I) 
ZI=ZO(I) 

IF(I.GE.NP2)  GOTO  788 
IF(LCHK.EQ.l)  GOTO  410 
JBEGIN=NABORS(I) 
JEND=NABORS (I+l) -1 
GOTO  415 
410       NABORS(I)  =K+1 
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JBEGIN=I+1 
JEND=NP 
415     CONTINUE 
C 
C   DECIDE  WHICH  MOLECULE  ATOM  I  IS  ON 

IM=(I-1)/NAM+1. 0  001 
C 
C   INNER  LOOP 

DO  290  JX=JBEGIN,JEND 
J=JX 

IF(LCHK.EQ.2)  J=LIST(JX) 
C 

C   IF  I  &  J  ARE  HEAD  GROUPS,  CALCULATE  HEAD-HEAD  REPULSIONS 
JCHK=J+NAM-1 

IF(MOD(ICHK,NAM) .NE. 0 . OR.MOD( JCHK,NAM) .NE.0)GOTO  614 
X=XI-XO(J) 
Y=YI-YO(J) 
Z=ZI-ZO(J) 
RSQ=X*X+Y*Y+Z*Z 
RIJ=SQRT(RSQ) 
RI=1.0/RIJ 
R3=RI*RI*RI 
R4=R3*RI 
R8=R4*R4 
C 
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ELJ=R3* (R8*RI+1 . 0) + (RIJ*HSHFTA) -HSHFT 
FLJS=3 . 0*R4* (4 . 0*R8*RI+1. 0) -HFSHFT 
FXD=FLJS*X*RI 
FYD=FLJS*Y*RI 
FZD=FLJS*Z*RI 
FX(I)=FX(I)+FXD 
FX(J)=FX(J)-FXD 
FY(I)=FY(I)+FYD 
FY(J)=FY(J)-FYD 
FZ(I)=FZ(I)+FZD 
FZ(J)=FZ(J)-FZD 
GO  TO  290 
C 
C   DECIDE  WHICH  MOLECULE  ATOM  J  IS  ON 

614    JM=(J-1)/NAM+1.0001 
C 
C   CHECK  WHETHER  I  AND  J  ARE  ON  SAME  MOL 

IF(IM.NE.JM)  GOTO  289 
C 

C   CHECK  IF  I  AND  J  ARE  FOURTH  NEIGHBORS  OR  GREATER 
KIJ=J-I 

IF(KIJ.LE.3)  GOTO  290 
C 

289   X=XI-XO(J) 
Y=YI-YO(J) 
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Z=ZI-ZO(J) 
C 

RSQ=X*X+Y*Y+Z*Z 
RIJ=SQRT(RSQ) 
IF(LCHK.NE.l)  GOTO  419 
IF(IM.EQ.JM)  GOTO  701 
IF(RSQ.GT.RDMAX)  GOTO  290 
IF(RSQ.GT.RLIST)  GOTO  290 
701      K=K+1 

LIST(K)=J 
419     IF(RIJ.GT.RC)  GOTO  290 

RI=1./RIJ 

R3=RI**3 

R6=R3*R3 

ELJ=2 . 0*R6* (R3-1 . 50) +ESHFT+RIJ*ESHFTA 

TOTLJ=TOTLJ+ELJ 

FLJS=18 . *RI*R6* (R3-1 . 0) -FSHFT 

FXD=FLJS*X*RI 

FYD=FLJS*Y*RI 

FZD=FLJS*Z*RI 

FX(I)=FX(I)+FXD 

FX(J)=FX(J)-FXD 

FY(I)=FY(I)+FYD 

FY(J)=FY(J)-FYD 

FZ(I)=FZ(I)+FZD 
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FZ(J)=FZ(J)-FZD 
290   CONTINUE 
C 

C   CALCULATE  BOND  LENGTHS 
C   CHECK  IF  LAST  ATOM  ON  MOLECULE 
788   IF(MOD(I,NAM) .EQ.O)  GOTO  300 
BX(I)=XI-X0(I+1) 
BY(I)=YI-Y0(I+1) 
BZ(I)=ZI-Z0(I+1) 

BSQ=BX(I) **2+BY(I) **2+BZ(I) **2 
B(I)=SQRT(BSQ) 
3  00     CONTINUE 
C 

IF(LCHK.NE.l)  GOTO  499 
NAB0RS(NP2)=K+1 
NABTOT=K 
C 
C   LOOP  CALCULATES  INTRAMOLECULAR  FORCES 

499   DO  600  1=1, NPl 
C 
C   CHECK  IF  LAST  ATOM  ON  MOLECULE 

IF (MOD (I, NAM) .EQ.O)  GOTO  600 
C 

C    CALCULATE  BOND  FORCES 
BDIFF=B(I)-BO 
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EBON=EBON+0 . 5*GAB*BDIFF*BDIFF 

FBON=GAB*BDIFF 

BI=1./B(I) 

BXI=BX(I)*BI 

BYI=BY(I) *BI 

BZI=BZ(I)*BI 

FXB=FBON*BXI 

FYB=FBON*BYI 

FZB=FBON*BZI 

FX(I)   =FX(I)-FXB 

FX(I+1)=FX(I+1)+FXB 

FY(I)   =FY(I)-FYB 

FY ( I+l ) =FY ( I+l ) +FYB 

FZ(I)   =FZ(I)-FZB 

FZ (I+l) =FZ (I+l) +FZB 
C 

C   CALCULATION  BENDING  FORCE 
C   CHECK  IF  NEXT-TO-LAST  ATOM  ON  MOLECULE 

IF(M0D(I+1,NAM) .EQ.O)  GOTO  600 

BII=1./B(I+1) 

BXJ=BX(I+1)*BII 

BYJ=BY(I+1)*BII 

BZJ=BZ(I+1)*BII 

CTI=- (BXJ*BXI+BYJ*BYI+BZJ*BZI) 

CDIFF=CTO-CTI 
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EBEN=EBEN+0 . 5*GTHA*CDIFF*CDIFF 

FT=GTHA*CDIFF 

FT1=FT*BI 

FT2=FT*BII 

FX1=FT1* (BXJ+CTI*BXI) 

FY1=FT1* (BYJ+CTI*BYI) 

FZ1=FT1* (BZJ+CTI*BZI) 

FX2=FT2* (BXI+CTI*BXJ) 

FY2=FT2* (BYI+CTI*BYJ) 

FZ2=FT2* (BZI+CTI*BZJ) 

FX(I)   =FX(I)   -FXl 

FX(H-1)=FX(I+1)+FX1-FX2 

FX ( 1+2 ) =FX ( 1+2 ) +FX2 

FY(I)   =FY(I)   -FYl 

FY(I+1)=FY(I+1)+FY1-FY2 

FY(I+2)=FY(I+2)+FY2 

FZ(I)   =FZ(I)   -FZl 

FZ (I+l) =FZ (I+l) +FZ1-FZ2 

FZ(I+2)=FZ(I+2)+FZ2 
C 

C    CALCULATION  OF  TORSION  FORCES 
C    CHECK  WHETHER  FIRST  ATOM  OF  MOLECULE 

IF(M0D(I+NAM-1,NAM) .EQ.O)  GOTO  600 

CALL  RYF0R(BX(I-1) ,BY(I-1) ,BZ(I-1) , BX (I) , BY (I) , BZ (I) , 
&       BX(I+1) ,BY(I+1) ,BZ(I+1) ,FX1,FY1,FZ1,FX2,FY2,FZ2, 
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&       FX3 , FYS , FZ3 , FX4 , FY4 , FZ4 , ETQ, CODI , GTHE) 
C 

IC0UNT=IC0UNT+1 

IF (CODI. LT. -0.50)  ITRAN=ITRAN+1 

FX(I-1)=FX(I-1)+FX1 

FY(I-1)=FY(I-1)+FY1 

FZ(I-1)=FZ(I-1)+FZ1 

FX(I)   =FX(I)   +FX2 

FY (I)   =FY(I)   +FY2 

FZ(I)   =FZ(I)   +FZ2 

FX(I+1)=FX(I+1)+FX3 

FY(I+1)=FY(I+1)+FY3 

FZ(I+1)=FZ(I+1)+FZ3 

FX ( 1+2 ) =FX ( 1+2 ) +FX4 

FY(I+2)=FY(I+2)+FY4 

FZ ( 1+2 ) =FZ ( 1+2 ) +FZ4 

ETOR=ETOR+ETQ 
600     CONTINUE 


C 


CALCULATE  GROUP-SHELL  INTERACTION 
EWALL=0.0 

RCHKSQ=(RWALL-2.0) **2 
DO  650  1=1, NP 
ICHK=I-1+NAM 
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XI=XO(I) 
YI=YO(I) 
ZI=ZO(I) 

c 

C   DISTANCE  FROM  CENTER  TO  GROUP 

RSQ=XI*XI+YI*YI+ZI*ZI 
C   IF  HEAD  GROUP,  CALCULATE  INTERACTIONS  WITH  WALL  USING  A 
C        POTENTIAL  OF  THE  FORM  —   GAMMA* (R-RO) **2 
IF(MOD(ICHK,NAM) .NE.O)  GO  TO  384 

RR=SQRT(RSQ) 

RI=1.0/RR 

RD=RWALL-RR 
C 

EWALL=EWALL+FORCON* (RD-REQUIL) * (RD-REQUIL) 

FWALL=-2 . 0*FORCON* (RD-REQUIL) 

FTOTWL=FTOTWL+FWALL 
C 

FX ( I ) =FX ( I ) -FWALL*XI *RI 

FY ( I ) =FY ( I ) -FWALL* YI *RI 

FZ (I) =FZ (I) -FWALL*ZI*RI 

GO  TO  650 
C  IF  NOT  A  HEAD  GROUP,  CALCULATE  REPULSIVE  INTERACTION  WITH 
C     THE  WALL 

384      IF(RSQ.LT.RCHKSQ)  GOTO  650 
C 


210 

C   DISTANCE  FROM  GROUP  TO  SHELL 

RR=SQRT(RSQ) 

RI=1.0/RR 

RD=RWALL-RR 

RDR=1.0/RD 

R3  =RDR*RDR*RDR 

R6  =R3*R3 

R12=R6*R6 
C 
C   FORCE  AND  ENERGY  FOR  R-12  WALL 

EWALL=EWALL+R12 

FWALL=12 . 0*R12*RDR 

FTOTWL=  FTOTWL+  FWALL 

FX ( I ) =FX ( I ) -FWALL*XI *RI 

FY ( I ) =FY ( I ) -FWALL*YI*RI 

FZ (I) =FZ (I) -FWALL*ZI*RI 
650    CONTINUE 
C 
C   CALCULATE  WALL  PRESSURE  IN  ATMOSPHERES 

EXTPRS=0.0085382*FTOTWL/RWALL**2 
C 

C   CALCULATE  INTRAMOLECULAR  AND  TOTAL  ENERGY 
EINTR=ETOR  +  EBON  +  EBEN 
TOTE  =  TOTLJ  +  EINTR  +  EWALL 
TRFRAC=FLOAT ( ITRAN) /FLOAT ( ICOUNT) 
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RETURN 
END 


SUBROUTINE  RYFOR ( BONXl , BONYl , BONZ 1 , B0NX2 , B0NY2 , B0NZ2 , 
&    B0NX3 , B0NY3 , B0NZ3 , FORXl , FORYl , FORZ 1 , F0RX2 , F0RY2 , 
&    FORZ 2 , F0RX3 , F0RY3 , FORZ 3 , F0RX4 , F0RY4 , FORZ 4 , ETOR, 
&    CODI,GR) 
C 

DATA  A1,A2,A3,A4,A5/-1.4  62,-1.578,0.3  68, 
&       3.156,3.788/ 
C 
C   CALCULATE  NORMAL  TO  PLANE  OF  BONDS  1  AND  2 

E=BONYl*BONZ2-BONZl*BONY2 

D=BONZl*BONX2-BONXl*BONZ2 

C-BONXl*BONY2-BONYl*BONX2 

AN0RM2=C*C+D*D+E*E 

ANORM  =SQRT(AN0RM2) 

AN0RMR=1 . /ANORM 

AX=E*ANORMR 

AY=D*ANORMR 

AZ=C*ANORMR 
C 
C   CALCULATE  NORMAL  TO  PLANE  OF  BONDS  2  AND  3 

H=BONY2*BONZ3-BONZ2*BONY3 
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G=BONZ  2  *B0NX3 -B0NX2  *BONZ  3 

F=B0NX2  *B0NY3 -B0NY2  *B0NX3 

BN0RM2=F*F+G*G+H*H 

BNORM  =SQRT(BN0RM2) 

BN0RMR=1 , /BNORM 

BX=H*BNORMR 

BY=G*BNORMR 

BZ=F*BNORMR 
C 
C   DETERMINE  COSINE  OF  THE  DIHEDRAL  ANGLE,  PHI 

CODI=  AX*BX+AY*BY+AZ*BZ 
C 
C   DETERMINE  TORSIONAL  ENERGY 

CODISQ=CODI*CODI 

C0DI3  =CODI*CODISQ 

ETOR   =  GR*(1.116+Al*CODI+A2*CODISQ+A3*CODI3 
&  +A4*CODISQ*CODISQ+A5*CODISQ*CODI3) 

C 
C   CALCULATE  DU/ COS (PHI) 

DUDCO= (Al+2 . *A2*CODI+3 . *A3*CODISQ+4 . *A4*C0DISQ*C0DI 
&         +5.*A5*C0DISQ*C0DISQ)*GR 
C 

C   START  EVALUATION  OF  FORCE  ON  GROUP  1 
C     CALCULATE  DANORM/DRl 

FACTl==(C*BONY2-D*BONZ2)/ANORM2 
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DADX1X=  -FACT1*AX 

DADX1Y=  -FACT1*AY-B0NZ2*AN0RMR 

DADX1Z=  -FACT1*AZ+B0NY2*AN0RMR 

C 

FACT2  =(E*BONZ2-C*BONX2)/ANORM2 
DADY1X=  -FACT2*AX+BONZ2*ANORMR 
DADY1Y=  -FACT2*AY 
DADY1Z=  -FACT2*AZ-BONX2*ANORMR 

C 

FACTS  =(D*BONX2-E*BONY2)/ANORM2 
DADZ1X=  -FACT3*AX-BONY2*ANORMR 
DADZ1Y=  -FACT3*AY+BONX2*ANORMR 
DADZ1Z=  -FACT3*AZ 

C 

C   CALCULATE  D(AB)/DR1 

DABDX1=BX*DADX1X+BY*DADX1Y+BZ*DADX1Z 
DABDY1=BX*DADY1X+BY*DADY1Y+BZ*DADY1Z 
DABDZ1=BX*DADZ1X+BY*DADZ1Y+BZ*DADZ1Z 

C 

C   CALCULATE  FORCE  ON  GROUP  1 
F0RX1=  -DUDC0*DABDX1 
F0RY1=  -DUDC0*DABDY1 
F0RZ1=  -DUDC0*DABDZ1 

C 

C   START  EVALUATION  OF  FORCE  ON  GROUP  2 
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C     CALCULATE  DAN0RM/DR2 

FACT4  =(D*B0NZ1-C*B0NY1)/AN0RM2 

DADX2X=  -FACT4*AX-DADX1X 

DADX2Y=  -FACT4*AY+B0NZ1*AN0RMR-DADX1Y 

DADX2Z=  -FACT4*AZ-B0NY1*AN0RMR-DADX1Z 
C 

FACT5  =(C*B0NX1-E*B0NZ1)/AN0RM2 

DADY2X=  -FACT5*AX-B0NZ1*AN0RMR-DADY1X 

DADY2Y=  -FACT5*AY-DADY1Y 

DADY2Z=  -FACTS *AZ+B0NX1*AN0RMR-DADY1Z 
C 

FACT 6  =(E*B0NY1-D*B0NX1)/AN0RM2 

DADZ2X=  -FACT6*AX+B0NY1*AN0RMR-DADZ1X 

DADZ2Y=  -FACT6*AY-B0NX1*AN0RMR-DADZ1Y 

DADZ2Z=  -FACT6*AZ-DADZ1Z 
C 
C   EVALUATE  D(BN0RM)/DR2 

FACT7  =(F*BONY3-G*BONZ3)/BNORM2 

DBDX2X=  -FACT7*BX 

DBDX2Y=  -FACT7*BY-BONZ3*BNORMR 

DBDX2Z=  -FACT7*BZ+BONY3*BNORMR 
C 

FACTS  =(H*BONZ3-F*BONX3)/BNORM2 
DBDY2X=  -FACT8*BX+BONZ3*BNORMR 
DBDY2Y=  -FACTS *BY 
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DBDY2Z=  -FACT8*BZ-BONX3*BNORMR 
C 

FACT9  =(G*BONX3-H*BONY3)/BNORM2 

DBDZ2X=  -FACT9*BX-BONY3*BNORMR 

DBDZ2Y=  -FACT9*BY+BONX3*BNORMR 

DBDZ2Z=  -FACT9*BZ 
C 
C   CALCULATE  D  C0S(PHI)/DR2 

DABDX2=BX*DADX2X+BY*DADX2Y+BZ*DADX2Z 
&         +AX*DBDX2X+AY*DBDX2Y+AZ*DBDX2Z 

DABDY2=BX*DADY2X+BY*DADY2Y+BZ*DADY2Z 
&         +AX*DBDY2X+AY*DBDY2Y+AZ*DBDY2Z 

DABDZ2=BX*DADZ2X+BY*DADZ2Y+BZ*DADZ2Z 
&         +AX*DBDZ2X+AY*DBDZ2Y+AZ*DBDZ2Z 
C 
C   FORCE  ON  GROUP  2 

F0RX2  =  -DUDC0*DABDX2 

F0RY2  =  -DUDC0*DABDY2 

F0RZ2  =  -DUDC0*DABDZ2 
C 

C   START  EVALUATION  OF  FORCE  ON  GROUP  4 
C     CALCULATE  D(BN0RM)/DR4 

FACTl  =(F*BONY2-G*BONZ2)/BNORM2 

DBDX4X=  -FACT1*BX 

DBDX4Y=  -FACT1*BY-B0NZ2*BN0RMR 
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DBDX4Z=  -FACT1*BZ+B0NY2*BN0RMR 
C 

FACT2  =(H*BONZ2-F*BONX2)/BNORM2 

DBDY4X=  -FACT2*BX+BONZ2*BNORMR 

DBDY4Y=  -FACT2*BY 

DBDY4Z=  -FACT2*BZ-BONX2*BNORMR 
C 

FACT3  =(G*BONX2-H*BONY2)/BNORM2 

DBDZ4X=  -FACTS *BX-B0NY2*BN0RMR 

DBDZ4Y=  -FACT3*BY+BONX2*BNORMR 

DBDZ4Z=  -FACT3*BZ 
C 
C   EVALUATE  D  C0S(PHI)/DR4 

DABDX4=AX*DBDX4X+AY*DBDX4Y+AZ*DBDX4Z 

DABDY4=AX*DBDY4X+AY*DBDY4Y+AZ*DBDY4Z 

DABDZ4=AX*DBDZ4X+AY*DBDZ4Y+AZ*DBDZ4Z 
C 
C   EVALUATE  FORCE  ON  GROUP  4 

F0RX4=  -DUDC0*DABDX4 

F0RY4=  -DUDC0*DABDY4 

F0RZ4=  -DUDC0*DABDZ4 
C 
C   EVALUATE  FORCE  ON  GROUP  3  BY  NEWTON'S  THIRD  LAW 

F0RX3  =  -FORX1-FORX2-FORX4 

F0RY3  =  -FORY1-FORY2-FORY4 
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F0RZ3  =  -FORZ1-FORZ2-FORZ4 

RETURN 
END 


BLOCK  DATA 

C0MM0N/DER/X2 (216),Y2(216),Z2(216),X3(216),Y3(216), 
&   Z3 (216) ,X4(216) ,Y4 (216) ,Z4(216) ,X5(216) ,Y5(216) , 
&   Z5(216) 

COMMON/FOR/ FX (216) , FY (216) , FZ (216) 

COMMON/DISP/DAX(216) ,DAY(216) ,DAZ(216) ,X0L(216) , 
&   Y0L(216) ,Z0L(216) 

COMMON/ PRO P/SUME , XSUM 

COMMON/MASS/MASS (216) 

DATA  X3,Y3,Z3,X4,Y4,Z4,X5,Y5,Z5/1944*0./ 

DATA  DAX, DAY , DAZ , FX, FY, FZ/1296*0 ./ 

DATA  MASS/2 16*1/ 

END 


APPENDIX  E 
ESTIMATION  OF  RUN  3  HEAD  GROUP  PARAMETERS 

In  Run  3  of  the  molecular  dynamics  simulations  of 
micelles,  the  head  group  was  given  the  mass  and  diameter  of 
a  sulfate  ion.   While  it  was  desirable  to  model  this  ion  as 
the  polar  head  due  to  its  common  usage  in  surfactants, 
information  on  its  intramolecular  potentials  is  incomplete. 
Data  pertaining  to  other  molecular  groups  was  used  to 
estimate  intramolecular  potential  parameters  necessary  to 
complete  the  model  of  the  sulfate  head  group  (O'Connell, 
1987)  . 

Muller  and  coworkers  have  measured  the  mean  bond 
lengths,  bond  angles,  and  the  vibrational  and  bending 
amplitudes  for  the  HS04~  ion  (Muller  and  Nagarajan,  1967; 
Muller  et  al . ,  1968).   The  values  are  assigned  to  the  bonds 
and  angles  in  Figure  E-la.   Blukis  et  al.  (1963)  have 
measured  the  C-0  bond  length  and  C-C-0  angle  of  dimethyl 
ether.   These  values  are  given  in  Figure  E-lb.   These  are 
used  in  the  estimation  of  bond  lengths,  angles  and 
amplitudes  for  the  0-0-8-03"  arrangement  of  the  sulfate  head 
group  and  the  alpha  carbon  of  the  surfactant  chain. 
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Figure  E-1.   a.)  Schematic  representation  of  the  the  HS04~ 
ion,  showing  mean  bond  lengths,  bond  angles,  and  vibrational 
and  bending  amplitudes  (Muller  and  Nagarajan,  1967;  Muller 
et  al.,  1968).   b.)  The  C-C-0  bonds  of  dimethyl  ether, 
showing  the  mean  angle  and  the  mean  length  of  the  C-0  bond 
(Blukis  et  al. ,  1963)  . 
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In  the  simulation  model,  the  head  group  is  represented 
as  a  homogeneous  sphere  connected  to  the  surfactant  chain 
through  a  single  bond.   The  center  of  mass  of  the  atoms  of 
the  sulfate  group  determines  the  length  and  angle  of  this 
bond.   In  Figure  E-2,  the  sulfate  group  and  alpha  carbon  are 
schematically  represented.   The  appropriate  bond  lengths  and 
angles  have  been  used  from  Figure  E-1.   The  center  of  mass 
of  the  sulfate  group  is  shown  in  Figure  E-2,  and  the  dashed 
line  connecting  it  to  the  alpha  carbon  is  the  bond  for  the 
model  head  group.   Its  length,  2.6  Angstrom,  and  angle,  146 
degrees,  are  indicated  in  the  figure. 

The  force  constant  for  the  harmonic  bond  vibration 
potential  used  in  the  simulation  is  determined  from  the  mean 
amplitude  of  the  vibration: 

y=  1  .68X  10"^-^  (E.l) 

a 

where  y  is  the  force  constant  in  mdyne/ Angstrom. 

^  is  the  reduced  mass  of  the  two  molecular  groups. 

a  is  the  mean  amplitude. 
For  the  head  group  bond  (dashed  line  in  Figure  E-2)  a  mean 
amplitude  of  vibration  is  calculated  from  its  components  to 
be  .078  Angstrom.   In  equation  E.l  this  yields  vibrational 
force  constant  for  the  head  group  bond  to  the  alpha  carbon 
of  .37  mdyne/ Angstrom. 
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Figure  E-2 .   The  sulfate  head  group,  shown  with  two  segments 
of  the  surfactant  chain.   Its  center  of  mass  is  designated 
by  the  "*"  and  the  dashed  line  connecting  it  to  the  alpha 
carbon  designates  the  bond  length  and  angle  of  the  head 
group  of  Run  3.   In  this  simulation,  the  sulfate  group  is 
modeled  by  a  single  sphere,  centered  at  the  sulfate  center 
of  mass,  with  intramolecular  potentials  derived  from  the 
sulfate  group. 
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The  force  constant  for  the  bond  angle  bending  potential 
is  proportional  to  the  mass  of  the  head  group.   It  is 
therefore  seven  times  greater  in  Run  3  than  in  the  other 
runs.   Note  that  although  Run  2  used  the  same  head  group 
mass,  only  the  mass  effect  was  investigated  in  that 
simulation;  all  intramolecular  potentials  were  the  same  as 
in  Run  1 . 

Rotational  potential  barriers  were  obtained  for 

molecular  groups  of  geometry  similar  to  that  of  the  sulfate 

group  (Cahill  et  al.,  1968).   They  are,  in  kcal/mole: 

methyl  formate  1.2 

propionaldehyde  2 . 3 

methyl  vinyl  ether  3.8 

cis-1-butene  4.0 

By  comparison  to  these  values  the  rotational  barrier  for  the 

larger,  more  polar  sulfate  group  is  estimated  to  be  5.0 

kcal/mole. 


APPENDIX  F 

MOLECULAR  DYNAMICS  PROGRAM  LISTING 

FOR  SIMULATION  3 


C  MOLECULAR  DYNAMICS  SOURCE  PROGRAM  FOR  A  SIMPLE 

C  MODEL  OF  A  MICELLE  ENCLOSED  IN  A  SHELL.  A 

C  FIFTH-ORDER  PREDICTOR-CORRECTOR  ALGORITHM  IS 

C  USED  TO  SOLVE  THE  EQUATIONS  OF  MOTION.  A 

C  SKELETAL  MODEL  DUE  TO  WEBER  IS  USED  FOR  EACH 

C  ALKANE  MOLECULE. 

C 

C  HEAD  GROUP  MASS  CAN  BE  SET.   SHIFTED-FORCE 

C  POTENTIALS  AND  NEIGHBOR  LIST  ARE  USED. 


C 


PARAMETER (MHEAD=1 ) 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 

COMMON/ DER/X2 (216) , Y2 (216) , Z2 (216) , X3 (216) , Y3 (216) , 
&   Z3(216) ,X4(216) ,Y4(216) ,Z4(216) ,X5(216) ,Y5(216) , 
&   Z5(216) 

COMMON/FOR/FX(216) ,FY(216) ,FZ(216) 

COMMON/DISP/DAX(216) , DAY (216) ,DAZ(216) ,X0L(216) , 
&   Y0L(216) ,Z0L(216) 

COMMON/NABLST/LIST( 12  000) ,NABORS(216) 

COMMON/ CONT/ GAB , CTO , GTHA , GTHE , BO , B2 I , Al , AO , EPS 

COMMON/ PARM/FO 2 , F12 , F32 , F42 , F52 

COMMON/ PROP/ SUME , XSUM 

COMMON/ PRO PA/RC , RLIST , FSHFT , ESHFT , ESHFTA , 
&   CUBE , CUBE2 , RDMAX , RCHH , HSHFT , HSHFTA , HFSHFT , REQUIL, 


224 


225 

COMMON/CONT/Al , AO , EPS 
COMMON/ FARM/ FO 2 , F12 , F32 , F42 , F52 
COMMON/ PRO P/SUME , XSUM 

C0MM0N/PR0PA/REQUIL,F0RC0N,HDSIG,SIG12 

COMMON/ENERGY/TOTE , TOTLJ , ETOR , EBON , EBEN , EINTR , TRFRAC , 
&   EXTPRS , FTOTWL 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB 
COMMON/SCAL/RSPHER,RWALL, DELRS, ISCALE 
COMMON/MASS/MASS (216) 

COMMON/ I NTRA/BO (216) ,GAB(216) ,CT0(216) ,GTHA(216) , 
&   GTHE(216) 
C 

C    FLAG,  IF  IFLAG.NE.l  USES  INTPOS  AND  RANDOM  VELOCITIES 
C      IF  IFLAG.EQ.l  XO,YO,ZO,AND  ALL  DERIVATIVES  ARE  READ 
C      FROM  PREVIOUS  RUN. 
C 

IFLAG=1 
C 

C  ===  SET  NUMBER  OF  PARTICLES  IN  PRIMARY  CELL 
NAM=9 
NM=24 
NP=NM*NAM 
PART=NP 
NP1=NP-1 
NP2=NP-2 
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NP22=.5*PART+.01 
C 

C  ===  SET  HEAD  GROUP  MASS,  BOND,  AND  FORCE  PARAMETERS 
DO  20  1=1, NP, NAM 
MASS(I)=MHEAD 
BO(I)=BLHEAD 
GAB ( I ) =GLHEAD 
CTO(I)=BAHEAD 
GTHA ( I ) =GAHEAD 
GTHE(I)=GTHEAD 
20     CONTINUE 
C 

C  ===  SET  RELATIVE  SIGMA  OF  HEAD-ALKANE  PAIR 
HDSIG=HDSGMA 
SIG12=(l.+HDSIG)/2 
C 

C  ===  SET  RADIUS  OF  SPHERE  ON  WHICH  HEAD  GROUPS  ARE  ATTACHED 
RSPHER=  6.000 
DELRS=0.0025 
ISCALE=0 

RWALL  =RSPHER+1.0 
REQUIL=1. 
FORCON=3  0. 
C 
C  ===  SET  VALUES  OF  PHYSICAL  CONSTANTS 
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XNAM=FLOAT (NAM) 

WTM0L=XNAM*14 . +2 . 

WTPART=WTMOL/ XNAM 

RSTAR=0.4 

EPS=419. 

AVO=6.0225E+23 

BOLZ=1.38054D-23 

EPSI=EPS/BOLZ/AVO 

THIRD=l./3. 

PI=3. 1415926535 
C 
C  ===  SET  DESIRED  FLUID  STATE  CONDITION 

TR=5.913 

T=TR*EPSI 

V0L=1. 333333 *PI*RSPHER**3 

DR=FLOAT (NM) /VOL 
C 
C  ===  SET  RUN  FLAGS  AND  PARAMETERS 

IFLG=-1 

KB=0 

KSAVE=10 

KSORT=10 

KWRITE=10 

MAXKB=4  0000 

XDIST=0.1 
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C 

C  ===  SET  TIME-STEP  AND  ITS  MULTIPLES 

DELTA=0. 00060 

DELSQ=DELTA* DELTA 

DELTSQ=.5*DELSQ 

TSTEP=SQRT(WTPART*RSTAR**2/EPS/1.E21) *1.D+12 

TT1=TSTEP*DELTA*1.D-12 
C 
C  ===  SET  PARAMETERS  IN  PREDICTOR-CORRECTOR  METHOD 

F02=3./16. 

F12=251./360. 

F32=ll./18. 

F42=l./6. 

F52=l./60. 
C 
C  ===  SCALE  FACTOR  FOR  VELOCITIES  DURING  EQUILIBRATION 

AHEAT=DELSQ*PART*3 . *TR 
C 
C  ===  WRITE  TAPE  HEADING 

OPEN(20,FILE='C872DAT2 ' , ACCESS= ' SEQUENTIAL' , 
&  FORM= ' UNFORMATTED ' ) 
WRITE (20)  NP,NAM,NM,DR,TR,EPS,RSTAR 
WRITE (20)  VOL, DELTA, TSTEP 
C 
C  ===  INITIALIZE  SUM  ACCUMMULATORS 
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XSUM=0 . 
SUME=0. 

c 

C  ===  PRINT  PARAMETERS 

WRITE(*,900) 
900    FORMAT (IHl///) 

WRITE(*,902) 
902    FORMAT(7X,49( '*•) ) 

WRITE(*,904) 
904    F0RMAT(7X, •*' ,T56, '*•) 

WRITE(*,906)  NM 
906    FORMAT ( 7 X, '*' ,2X, 'MOLECULAR  DYNAMICS  FOR', 13, 
A  '  N-ALKANE  MOLECULES ', T56 ,'*' ) 

WRITE (*, 904) 

WRITE (*, 917)  NAM 
917    F0RMAT(7X, '*' ,5X, 'WITH' ,13, '  PARTICLES  PER  MOLECULE ' 
&   ,T56, '*') 

WRITE (*, 904) 

WRITE (*, 902) 

WRITE (*, 904) 

WRITE (*, 908)  EPSI,RSTAR 
908    FORMAT ( 7X, '*' ,2X, 'EPSI/K  =  ' , F7 . 3 , T36 , 'RSTAR  =  ', 
&    F7.3,T56, '*•) 

WRITE(*,910)  TR,T 
910    F0RMAT(7X, '*' ,2X, 'TR      =  ',F7.3,T36,'     T=  ', 
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&    F7.3,T56, '*') 

WRITE(*,912)  DR,VOL 
912    FORMAT ( 7X, '*• ,2X, 'DR      =  ' , F7 . 3 , T36 , ' VOL    =  ', 
&    F8.3,T56, '*' ) 

WRITE (*, 92  0)  DELTA 
920    FORMAT ( 7X, '*' ,2X, 'DELTA   =  ' , F9 . 5 ,T56 , ' * ' ) 
WRITE (*, 904) 
WRITE (*, 922)  TSTEP 
922    F0RMAT(7X,  '*' ,2X,  'TIME  UNIT  =  •,F6.3,'E-12  SEC, 
&    T56, '*•) 

WRITE(*,924)  TTl 
924    F0RMAT(7X, '*• ,2X, 'TIME  STEP  =  ',1PE10.3,'  SEC, 
&    T56, '*' ) 

WRITE (*, 904) 
WRITE (*, 926)  DE 
926    F0RMAT(7X, '*' ,2X, 'ENERGY    CORRECTION  =  ' , F7 . 3 , 
&    T56, '*') 

WRITE (*, 904) 
WRITE (*, 902) 
WRITE (*, 93  2)  RSPHER 
932     FORMAT (////I OX, 'RADIUS  OF  SPHERE  =  ',F6.2) 
WRITE(*,934)  GTHE(2) 


934     FORMAT (//lOX, 'ROTATIONAL  POTENTIAL  PARAMETER 
&    F7.3/) 


I 


r 
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C 

OPEN (30, FILE= ' LASTDAT ' , ACCESS= ' SEQUENTIAL ' , 
&    FORM= ' UNFORMATTED ' ) 

IF(IFLAG.EQ.l)  GOTO  199 
C 

C  ===  LOAD  INITIAL  POSITIONS  OF  ATOMS 
CALL  INTPOS(RSPHER) 
CALL  POSPRI(KB,NM,NAM) 
C 

C  ===  PRINT  RUN-TABLE  HEADING 
WRITE(*,930) 
93  0    FORMAT (1H1////4X, 'KB' ,5X, 'RSPH' ,4X, ' ENRG ' ,5X, 'El' , 
A    4X, 'DIST' ,6X, 'TEMP' ,3X, 'TRFRC ,3X, 'TOT  ENR', 
A    3X, 'FTOTWL' ,3X, 'PRESSURE'/) 
C 
C  ===  LOAD  INITIAL  VELOCITIES  OF  ATOMS 

CALL  INTVEL(AHEAT,PART) 
C 

C  ===  ASSIGN  INITIAL  ACCELERATIONS  BASED  ON  INITIAL 
C       POSITIONS 

CALL  EVAL(RWALL) 
C 
C  ===  SCALE  ACCELERATIONS  AND  STORE  STARTING  POSITIONS 
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DO  530  1=1, NP 

X2 ( I ) =FX ( I ) *  DELTSQ/MASS ( I ) 
Y2 ( I ) =FY ( I ) *DELTSQ/MASS ( I ) 
Z2 ( I ) =FZ ( I ) *DELTSQ/MASS ( I ) 
53  0    CONTINUE 
GOTO  188 
C 

C   READ  POSITIONS  AND  DERIVATIVES  INSTEAD  OF  USING 

C   FCC,  INVEL,  ACCELERATION  

C 

199  WRITE(*,930) 

READ(30)  KB,IFLG,NM,NAM,RSPHER,XSUM,SUME 
RWALL=RSPHER+1 . 0 
DO  200  K=1,NP 

READ(30)  X0(K) ,Y0(K) ,ZO(K) 

READ(30)  X1(K) ,Y1(K) ,Z1(K) 

READ(30)  X2(K) ,Y2 (K) ,Z2 (K) 

READ(30)  X3(K) ,Y3(K) ,Z3 (K) 

READ(30)  X4(K) ,Y4(K) ,Z4(K) 

READ(30)  X5(K) ,Y5(K) ,Z5(K) 

READ(30)  DAX(K) ,DAY(K) ,DAZ(K) 

200  CONTINUE 
REWIND  30 

188      DO  377  1=1, NP 
X0L(I)=X0(I) 
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YOL(I)=YO(I) 

ZOL(I)=ZO(I) 
377        CONTINUE 
C 

c 

c 

C  ===  ENTER  MAIN  LOOP  OF  SIMULATION 

IF(IFLAG.NE.l)  GOTO  777 

CALL  PREDCT(NP) 

CALL  EVAL(RWALL) 

CALL  CORR(DELTSQ) 

IF ( ISCALE . EQ . 1)  RSPHER=RSPHER+DELRS 

RWALL=RSPHER+1 . 0 
777       NS=KB+1 

DO  599  NTIMES=NS , MAXKB 

KB=KB+1 

CALL  PREDCT(NP) 

CALL  EVAL(RWALL) 

CALL  CORR(DELTSQ) 
C 
C  ===  CALCULATE  MEAN  SQUARE  DISPLACEMENT  &  KINETIC  ENERGY 

SUMVEL=0 . 

TDIST=0. 

DO  540  1=1, NP 

TDIST=TDIST+DAX(I)**2+DAY(I)**2+DAZ(I)**2 
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SUMVEL=SUMVEL+(X1(I)**2+Y1(I)**2+Z1(I)**2) * 
&        MASS (I) 
54  0      CONTINUE 

TDIST=TDIST/PART 
EK=SUMVEL/ (2 . *PART*DELSQ) 
C 

C  ===  ACCUMMULATE  SUMS  FOR  PROPERTY  AVERAGES 
XSUM=XSUM+SUMVEL 
SUME=SUME+TOTE 
C 

C  ===  PROPERTY  CALCULATION  &  PRINT-OUT  AT  INTERVALS 
IF(MOD(KB,KWRITE) .NE.O)  GOTO  550 
FKB=FLOAT (KB) *PART 
TMP=XSUM/ ( 3 . *DELSQ*FKB) 
ENR=(SUME/FKB) 
E1=T0TE/PART 
ET0T=E1+EK 

RLTIM= DELTA* FLOAT (KB) *TSTEP 

WRITE ( * , 9  4  0 )  KB , RSPHER , ENR , El , TDIST , TMP , TRFRAC , 
&        ETOT,FTOTWL,EXTPRS 
C 

IF(MOD(KB, 1000) .NE.O)  GOTO  550 
CALL  POSPRI(KB,NM,NAM) 
WRITE(*,930) 
940        F0RMAT(1H  , 16 , F6 . 3 , 3F9 . 2 , FIO . 3 , IX, F6 . 3 , IX, F9 . 3 , 
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&  2(2X,F9.3)) 

C 

C  ===  DURING  FIRST  OF  RUN,  SCALE  VELOCITIES  FOR  TEMPERATURE 
550      IF(IFLG.LT.l)  CALL  EQBRAT(SUMVEL,AHEAT,TDIST, 
&      XDIST,NP,IFLG,KB,NAM) 
C 

C  WRITE  DATA  ONTO  TAPE  FOR  LATER  USE  

IF(KB.EQ.O)  GOTO  777 
C        IF(IFLG.EQ.O)  GOTO  588 

IF(MOD(KB,10) .NE.O)  GOTO  588 

WRITE (20)  KB, TMP, ETOT, El , EK, TOTE , EINTR, TOTLJ , ETOR, 
&    EBON,EBEN 

WRITE(20)  XO,YO,ZO 
WRITE(20)  X1,Y1,Z1 
WRITE(20)  X2,Y2,Z2 
C         WRITE(20)  X3,Y3,Z3 
C         WRITE (20)  X4,Y4,Z4 
C         WRITE(20)  X5,Y5,Z5 

588    IF(MOD(NTIMES, 100) .NE.O)  GOTO  599 

WRITE (30)  KB,IFLG,NM,NAM,RSPHER,XSUM,SUME 
DO  598  K=1,NP 

WRITE(30)  XO(K) ,YO(K) ,ZO(K) 
WRITE(30)  X1(K) ,Y1(K) ,Z1(K) 
WRITE(30)  X2(K) ,Y2(K) ,Z2(K) 
WRITE(30)  X3(K) ,Y3(K) ,Z3(K) 
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WRITE(30)  X4(K) ,Y4(K) ,Z4(K) 
WRITE(30)  X5(K) ,Y5(K) ,Z5(K) 
WRITE(30)  DAX(K) ,DAY(K) ,DAZ(K) 

598  CONTINUE 
REWIND  3  0 

599  CONTINUE 
CLOSE (20) 
CLOSE (3  0) 
STOP 

END 

SUBROUTINE  POSPRI  (KB,  N]yi,NAM) 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

11=0 

WRITE (*, 942)  KB 
942      FORMAT (1H1///7X, 'POSITIONS  OF  GROUPS  AT  TIME-STEP:' 
&     ,16//) 

DO  404  JJ=1,NM 
WRITE(*,935) 
93  5         FORMAT (//2X, 'MOLECULE' , 7X, ' X ' , IIX, ' Y ' ,12X, 'Z' , 
&       12X, 'R' ,9X, 'BOND  LEN'/) 
DO  4  04  KK=1,NAM 
11=11+1 

RR=SQRT(X0(II)**2+Y0(II)**2+Z0(II)**2) 
XXX=X0 (II) -XO (II-l) 
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YYY=Y0(II)-Y0(II-1) 
ZZZ=Z0(II)-Z0(II-1) 
BB=SQRT(XXX*XXX+YYY*YYY+ZZZ*ZZZ) 
IF(KK.EQ.l)  WRITE(*,909)  JJ, XO (II) , YO (II) , 
&  ZO(II) ,RR 

IF(KK.GT.l)  WRITE(*,911)  XO (II) , YO (II) , ZO (II) , 
&  RR,BB 

4  04      CONTINUE 

909      FORMAT(5X,I2,4 (5X,F8.4)  ) 
911      FORMAT(7X,5(5X,F8.4)  ) 
RETURN 
END 


SUBROUTINE  INTPOS (RSPHER) 
C 

C   SET-UP  INITIAL  POSITIONS  OF  PARTICLES.  HEAD  GROUPS  ARE 
C     ASSIGNED  TO  FIXED  POSITIONS  ON  THE  SURFACE  OF  A  SPHERE 
C     OF  RADIUS  "RSPHER". 


C 


COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 
COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB 
COMMON/CONT/Al , AO , EPS 

COMMON/INTRA/B0(216) , GAB (2 16) ,CT0(216) ,GTHA(216) , 
GTHE(216) 
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DIMENSION  THETA(50) ,PHI(50) 
C 

PI=3. 1415926535 
PI180=PI/180. 
C 

C   TWENTY-FOUR  MOLECULES 
THETA(1)=0. 
PHI(1)=0. 
THETA(24)=PI 
PHI(24)=0. 
PH24=0. 
DO  70  1=2,5 

PHI(I)=PH2  4 
PH24=PH24+PI/2. 

70  THETA(I)=36.*PI180 
DELPHI=2.*PI/7 
PH24=DELPHI/2. 

DO  71  1=5,12 
PHI(I)=PH24 
PH2  4=PH2  4+DELPHI 

71  THETA(I)=72.*PI180 
PH24=0. 

DO  72  1=13,19 
PHI(I)=PH2  4 
PH24=PH24+DELPHI 
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72  THETA(I)=108.*PI180 
PH24=PI/4. 

DO  73  1=20,23 
PHI(I)=PH2  4 
PH24=PH24+PI/2. 

73  THETA(I)=144.*PI180 
C 

60    WRITE(*,909) 
909     FORMAT (1H1////7X, 'INITIAL  POSITIONS  OF  HEAD  GROUPS' 
&    //7X, 'MOLECULE' ,5X, ' THETA ' ,8X, 'PHI'/) 
DO  101  1=1, NM 

TTT=THETA ( I ) /PI 180 
PPP=PHI(I)/PI180 
WRITE (*, 900)  I,TTT,PPP 
900       FORMAT(10X,I3,2(6X,F7.3) ) 
101     CONTINUE 
C 

c 

C    STANDARD  DISTANCES  ALONG  TRANS  CHAIN 

XLEN=0.5*B0(2) *SQRT (2 . * ( 1 . -CTO (2 ) ) ) 

ZLEN=SQRT ( BO ( 2 ) *B0 ( 2 ) -XLEN*XLEN) 

ZLENSQ=ZLEN*ZLEN 
C 
C    HEAD-TO-ALPHA  DISTANCE 

XLEN1=0.5*B0(1) *SQRT(2.*(1.-CT0(1) ) ) 
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ZLEN1=SQRT(B0(1)*B0(1)-XLEN1*XLEN1) 
ZLN1SQ=ZLEN1*ZLEN1 
C 

C  ===   ADD  ATOMS  TO  FORM  MOL  IN  TRANS  POS  = 
DO  120  I-1,NM 
IADD=(I-1)  *NA]yi 
THR=THETA ( I ) 
PHIR=PHI ( I ) 
SR=SIN(THR) 
CR=COS(THR) 
SPR=SIN(PHIR) 
CPR=COS(PHIR) 
C 

C   HEAD  GROUP  ON  A  MOLECULE 
J=l 

K=IADD+J 
RK=RSPHER 
XO(K)=RK*SR*CPR 
YO(K)=RK*SR*SPR 
ZO(K)=RK*CR 
C 

C   ALPHA  GROUP  ON  A  MOLECULE 
J=2 

K=IADD+J 
RDISA=RSPHER-XLEN1 
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RK=SQRT (RDISA*RDISA+ZLN1SQ) 

BETA=ASIN(ZLEN1/RK) 

THE=THR-BETA 

ST=SIN(THE) 

XO(K)=RK*ST*CPR 

YO(K)=RK*ST*SPR 

Z0(K)=RK*COS(THE) 
C 

C   ODD  NUMBERED  GROUPS  ON  A  MOLECULE 
DO  92  J=3,NAM,2 

K=IADD+J 

RK=RDISA-FLOAT ( J-2 ) *XLEN 

XO(K)=RK*SR*CPR 

YO(K)=RK*SR*SPR 

ZO(K)=RK*CR 
92       CONTINUE 
C 

C   EVEN  NUMBERED  GROUPS  ON  A  MOLECULE 
DO  9  5  J=4,NAM,2 

K=IADD+J 

RDIS=RDISA-FLOAT (J-2 ) *XLEN 

RK=SQRT (RDIS*RDIS+ZLENSQ) 

BETA=ASIN(ZLEN/RK) 

THE=THR-BETA 

ST=SIN(THE) 
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XO(K)=RK*ST*CPR 
YO(K)=RK*ST*SPR 
Z0(K)=RK*COS(THE) 
9  5       CONTINUE 
12  0     CONTINUE 
RETURN 
END 


SUBROUTINE  INTVEL ( AHEAT , PART) 
C 

C  ===  ASSIGN  INITIAL  VELOCITIES  TO  ATOMS 
C 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB 
COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 
COMMON/MASS/MASS (216) 
SUMX=0 . 
SUMY=0 . 
SUMZ=0. 
DO  200  1=1, NM 
XX=RANF(DUM) 
YY=RANF(DUM) 
ZZ=RANF(DUM) 
NX=(I-1)*NAM 
DO  2  00  J=1,NAM 
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X1(J+NX)=XX 
Y1(J+NX)=YY 
Z1(J+NX)=ZZ 

SUMX=SUMX+X1 (J+NX) *MASS ( J+NX) 
SUMY=SUMY+Y1 (J+NX) *MASS (J+NX) 
SUMZ=SUMZ+Z1 (J+NX) *MASS (J+NX) 
2  00    CONTINUE 
C 

C  ===  SCALE  VELOCITIES  SO  THAT  TOTAL  MOMENTUM  =  ZERO 
X=0. 
DO  210  1=1, NP 

XI ( I ) =X1 ( I ) -SUMX/PART/MASS ( I ) 
Yl ( I ) =Y1 ( I ) -SUMY/ PART/MASS ( I ) 
Z 1 ( I ) =Z 1 ( I ) -SUMZ/ PART/MASS ( I ) 
X=X+(X1(I)**2+Y1(I)**2+Z1(I)**2)*MASS(I) 
210    CONTINUE 
C 

C  ===  SCALE  VELOCITIES  TO  DESIRED  TEMPERATURE 
HEAT=SQRT (AHEAT/X) 
DO  220  1=1, NP 

X1(I)=X1(I)*HEAT 
Y1(I)=Y1(I)*HEAT 
Z1(I)=Z1(I) *HEAT 
220    CONTINUE 
RETURN 
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END 


C 

c  == 
c 
c 


SUBROUTINE  PREDCT(NP) 

=  USE  TAYLOR  SERIES  TO  PREDICT  POSITIONS  &  THEIR 
DERIVATIVES  AT  NEXT  TIME-STEP 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 

COMMON/DER/X2(216) , Y2 (216) , Z2 (216) , X3 (216) , Y3 (216) , 
&  Z3(216) ,X4(216) ,Y4(216) ,Z4(216) ,X5(216) ,Y5(216) , 
&     Z5(216) 

CO]yiMON/FOR/FX(216)  ,  FY  (216)  ,  FZ  (216) 
DO  300  1=1, NP 


XO(I 
YO(I 
ZO(I 
X1(I 
Y1(I 
Z1(I 
X2(I 
Y2(I 
Z2(I 
X3(I 


=X0 ( I ) +X1 ( I ) +X2 ( I ) +X3 ( I ) +X4 ( I ) +X5 ( I ) 

=Y0(I)+Y1(I)+Y2(I)+Y3(I)+Y4(I)+Y5(I) 

=Z0(I)+Z1(I)+Z2(I)+Z3(I)+Z4(I)+Z5(I) 

=X1(I)+2.*X2 (I)+3.*X3(I)+4.*X4 (I)+5.*X5(I) 

=Y1(I)+2.*Y2(I)+3.*Y3(I)+4.*Y4(I)+5.*Y5(I) 

=Z1(I)+2.*Z2(I)+3.*Z3(I)+4.*Z4(I)+5.*Z5(I) 

=X2(I)+3.*X3 (I)+6.*X4 (I)+10.*X5(I) 

=Y2(I)+3.*Y3(I)+6.*Y4(I)+10.*Y5(I) 

=Z2(I)+3.*Z3(I)+6.*Z4(I)+10.*Z5(I) 

=X3(I)+4.*X4(I)+10.*X5(I) 
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300 


Y3(I 

Z3(I 

X4(I 

Y4(I 

Z4(I 

FX(I 

FY  (I 

FZ(I 
CONTINUE 
RETURN 
END 


=Y3 (I)+4.*Y4 (I)+10.*Y5(I) 

=Z3(I)+4.*Z4(I)+10,*Z5(I) 

=X4(I)+5.*X5(I) 

=Y4(I)+5.*Y5(I) 

=Z4(I)+5.*Z5(I) 

=0. 

=  0. 

=  0. 


SUBROUTINE  CORR(DELTSQ) 
C 

C  ===  CORRECT  PREDICTED  POSITIONS  AND  THEIR  DERIVATIVES 
C 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 
COMMON/VEL/Xl(216) , Yl (216) , Zl (216) 

COMMON/DER/X2(216) , Y2 (216) , Z2 (216) , X3 (216) , Y3 (216) , 
&  Z3(216),X4(216),Y4(216),Z4(216),X5(216),Y5(216), 
&     Z5(216) 

COMMON/FOR/FX(216) , FY (216) , FZ (216) 

COMMON/DISP/DAX(216) ,DAY(216) ,DAZ(216) ,X0L(216) , 
&     Y0L(216) ,Z0L(216) 


246 

COMMON/ FARM/ FO 2 , F12 , F32 , F42 , F52 
COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB 
COMMON/PROP/SUME , XSUM 

COMMON/SCAL/RSPHER, RWALL, DELRS , ISCALE 
COMMON/MASS/MASS (216) 

IF(MOD(KB, 10) .NE. 0. OR. ISCALE. NE.l)  GOTO  610 
IF(RSPHER  .LE.  3.6)  GO  TO  610 
RS  PHER=RS  PHER- DELRS 
RWALL  =RWALL  -DELRS 


610    DO  690  1=1, NP 

XERR=X2 ( I ) -DELTSQ*FX ( I ) /MASS ( I ) 
YERR=Y2 ( I ) -DELTSQ*FY ( I ) /MASS ( I ) 
ZERR=Z2 ( I ) -DELTSQ*FZ ( I ) /MASS ( I ) 


X0(I 
X1(I 
X2(I 
X3(I 
X4(I 
X5(I 
Y0(I 
Y1(I 
Y2(I 
Y3(I 


=XO(I 
=X1(I 
=X2(I 
=X3(I 
=X4  (I 
=X5(I 
=YO(I 
=Y1(I 
=Y2(I 
=Y3(I 


-XERR*F02 

-XERR*F12 

-XERR 

-XERR*F32 

-XERR*F42 

-XERR*F52 

-YERR*F02 

-YERR*F12 

-YERR 

-YERR*F32 
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Y4(I)=Y4 (I)-YERR*F42 
Y5(I)=Y5(I) -YERR*F52 
Z0(I)=Z0(I)-ZERR*F02 
Z1(I)=Z1(I)-ZERR*F12 
Z2(I)=Z2(I)-ZERR 
Z3 (I)=Z3 (I) -ZERR*F32 
Z4 (I)=Z4(I)-ZERR*F4  2 
Z5(I)=Z5(I)-ZERR*F52 
690   CONTINUE 

IF(MOD(KB,KSORT) .NE.O.OR.ISCALE.NE.l)  GOTO  680 
DO  691  1=1, NP 

RSQ=XO(I) *XO(I)+YO(I)*YO(I)+ZO(I)*ZO(I) 
R=SQRT(RSQ) 
RFACT  =  1.0-DELRS/R 
XO(I)  =  XO(I)*RFACT 
YO(I)  =  YO(I)*RFACT 
ZO(I)  =  ZO(I)*RFACT 
691    CONTINUE 
C 

C  ===  DISPLACEMENTS 
680   DO  692  1=1, NP 

DAX ( I ) =DAX ( I ) -XO ( I ) +XOL ( I ) 
DAY(I)=DAY(I)-YO(I)+YOL(I) 
DAZ(I)=DAZ(I)-ZO(I)+ZOL(I) 
C 
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C  ===  STORE  NEW  POSITIONS 
XOL(I)=XO(I) 
YOL(I)=YO(I) 
ZOL(I)=ZO(I) 
692    CONTINUE 
RETURN 
END 


SUBROUTINE  EQBRAT ( SUMVEL, AHEAT , TDIST , XDIST , NP , IFLG , 
&   KB, NAM) 
C 

C  ===  SCALE  VELOCITIES  DURING  INITIAL  TIME-STEPS 
C 

COMMON/VEL/Xl(216) ,Y1(216) ,Z1(216) 
COMMON/ PROP/ SUME , XSUM 

COMMON/SCAL/RSPHER,RWALL,DELRS,ISCALE 
C 

IF(IFLG.EQ.-l)  GOTO  720 

IF(TDIST.GT.XDIST.0R.IFLG.EQ.-2)  GOTO  750 
72  0    HEAT=SQRT(AHEAT/SUMVEL) 
DO  730  1=1, NP 

X1(I)=X1(I) *HEAT 
Y1(I)=Y1(I)*HEAT 
Z1(I)=Z1(I)*HEAT 
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73  0    CONTINUE 

RETURN 
C 

C  ===  AT  END  OF  EQUILIBRATION  STAGE,  SET  PROPERTY  SUMS 
C       TO  ZERO 
750    IFLG=  1 

KB=0 

SUME=0. 

XSUM=0 . 

RETURN 

END 


SUBROUTINE  EVAL(RWALL) 

COMMON/POS/X0(216) ,Y0(216) ,Z0(216) 

COMMON/FOR/FX(216) ,FY(216) ,FZ(216) 

COMMON/ CONT/ A 1 , AO , EPS 

COMMON/ PROP/ SUME , XSUM 

C0MM0N/PR0PA/REQUIL,F0RC0N,HDSIG,SIG12 

COMMON/ENERGY/TOTE , TOTLJ , ETOR , EBON , EBEN , EINTR , 
&    TRFRAC , EXTPRS , FTOTWL 

COMMON/NUM/NM , NAM , NP , NPl , NP2 , NP2  2 , KSORT , KB 

COMMON/INTRA/B0(216) ,GAB(216) ,CT0(216) ,GTHA(216) , 
&    GTHE(216) 
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DIMENSION  B(216) ,BX(216) ,BY(216) ,BZ(216) 
C 

C   INITIALIZE  ACCUMULATORS 
ICOUNT=0 
K=0 

ITRAN  =0 
TOTE=0 . 0 
TOTLJ=0.0 
FTOTWL=0.0 
ETOR=0 . 0 
EBON=0 . 0 
EBEN=0.0 
EINTR=0.0 
DO  150  1=1, NP 
FX(I)=0.0 
FY (I) =0.0 
FZ(I)=0.0 
150   CONTINUE 
C 

C   CALCULATE  LJ  FORCES 
C 

C   OUTER  LOOP  OVER  PARTICLES 
DO  300  1=1, NPl 
ICHK=I-1+NAM 
XI=XO(I) 
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YI=YO(I) 
ZI=ZO(I) 

IF(I.GE.NP2)  GOTO  788 
C 
C   DECIDE  WHICH  MOLECULE  ATOM  I  IS  ON 

IM=(I-1)/NAM+1.0001 
C 
C   INNER  LOOP 

DO  290  J=I+1,NP 
C 

C   IF  I  &  J  ARE  HEAD  GROUPS,  CALCULATE  HEAD-HEAD  REPULSIONS 
JCHK=J+NAM-1 

IF(MOD(ICHK,NAM) . NE. 0 . OR.MOD ( JCHK, NAM) .NE.0)GOTO  614 
X=XI-XO(J) 
Y=YI-YO(J) 
Z=ZI-ZO(J) 
RSQ=X*X+Y*Y+Z*Z 
RIJ=SQRT(RSQ) 
RI=HDSIG/RIJ 
R3=RI*RI*RI 
R4=R3*RI 
R8=R4*R4 

HS3=1.0/HDSIG**3 
HS4=HS3/HDSIG 
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ELJ=R3* (R8*RI+HS3) 
FLJS=3 . 0*R4* (4 . 0*R8*RI+HS4) 
FXD=FLJS*X*RI 
FYD=FLJS*Y*RI 
FZD=FLJS*Z*RI 
FX(I)=FX(I)+FXD 
FX(J)=FX(J)-FXD 
FY(I)=FY(I)+FYD 
FY(J)=FY(J)-FYD 
FZ(I)=FZ(I)+FZD 
FZ(J)=FZ(J)-FZD 
GO  TO  290 
C 
C   DECIDE  WHICH  MOLECULE  ATOM  J  IS  ON 

614    JM=(J-1)/NAM+1.0001 
C 
C   CHECK  WHETHER  I  AND  J  ARE  ON  SAME  MOL 

IF(IM.NE.JM)  GOTO  289 
C 

C   CHECK  IF  I  AND  J  ARE  FOURTH  NEIGHBORS  OR  GREATER 
KIJ=J-I 

IF(KIJ.LE.3)  GOTO  290 
C 

289   X=XI-XO(J) 
Y=YI-YO(J) 
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Z=ZI-ZO(J) 
C 

RSQ=X*X+Y*Y+Z*Z 
RIJ=SQRT(RSQ) 
C 

C   IF  ONE  IS  A  HEAD  GROUP,  CALCULATE  HEAD-ALKANE  INTERACTION 
IF(MOD(ICHK,NAM) . EQ. 0 . OR.MOD ( JCHK, NAM) .EQ.O)  THEN 
RI=SIG12/RIJ 
R3=RI**3 
R6=R3*R3 

ELJ=2 . 0*R6* (R3-1. 50) 
FLJ=18 . *RI*R6* (R3-1 .0) 
C 

C   IF  NOT,  CALCULATE  ALKANE-ALKANE  INTERACTION 
ELSE 

RI=1./RIJ 

R3=RI**3 

R6=R3*R3 

ELJ=2 . 0*R6* (R3-1 .  50) 

FLJS=18 . *RI*R6* (R3-1 . 0) 
ENDIF 

TOTLJ=TOTLJ+ELJ 
FXD=FLJS*X*RI 
FYD=FLJS*Y*RI 
FZD=FLJS*Z*RI 
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FX(I)=FX(I)+FXD 

FX(J)=FX(J)-FXD 

FY(I)=FY(I)+FYD 

FY(J)=FY(J)-FYD 

FZ(I)=FZ(I)+FZD 

FZ(J)=FZ(J)-FZD 
290   CONTINUE 
C 

C   CALCULATE  BOND  LENGTHS 
C   CHECK  IF  LAST  ATOM  ON  MOLECULE 
788   IF(MOD(I,NAM) .EQ.O)  GOTO  300 

BX(I)=XI-X0(I+1) 

BY(I)=YI-Y0(I+1) 

BZ(I)=ZI-Z0(I+1) 

BSQ=BX(I)**2+BY(I) **2  +  BZ(I)  **2 

B(I)=SQRT(BSQ) 
300     CONTINUE 
C 

c 

C   LOOP  CALCULATES  INTRAMOLECULAR  FORCES 

499   DO  600  1=1, NPl 
C 
C   CHECK  IF  LAST  ATOM  ON  MOLECULE 

IF (MOD (I, NAM) .EQ.O)  GOTO  600 
C 
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C    CALCULATE  BOND  FORCES 

BDIFF=B(I)-BO(I) 

EBON=EBON+0.5*GAB(I) *BDIFF*BDIFF 

FBON=GAB ( I ) *BDIFF 

BI=1./B(I) 

BXI=BX(I)*BI 

BYI=BY(I)*BI 

BZI=BZ(I) *BI 

FXB=FBON*BXI 

FYB=FBON*BYI 

FZB=FBON*BZI 

FX(I)   =FX(I)-FXB 

FX(I+1)=FX(H-1)+FXB 

FY(I)   =FY(I)-FYB 

FY(I  +  1)=FY(I+1)+FYB 

FZ(I)   =FZ(I)-FZB 

FZ(I+1)=FZ(I+1)+FZB 
C 

C   CALCULATION  BENDING  FORCE 
C   CHECK  IF  NEXT-TO-LAST  ATOM  ON  MOLECULE 

IF(M0D(I+1,NAM) .EQ.O)  GOTO  600 

BII=1./B(I+1) 

BXJ=BX(I+1)*BII 

BYJ=BY(I+1) *BII 

BZJ=BZ(I+1)*BII 
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CTI=-  (BXJ*BXI+BYJ*BYH-BZJ*BZI) 

CDIFF=CTO(I)-CTI 

EBEN=EBEN+0.5*GTHA(I) *CDIFF*CDIFF 

FT=GTHA(I)*CDIFF 

FT1=FT*BI 

FT2=FT*BII 

FX1=FT1* (BXJ+CTI*BXI) 

FY1=FT1* (BYJ+CTI*BYI) 

FZ1=FT1* (BZJ+CTI*BZI) 

FX2=FT2*  (BXH-CTI*BXJ) 

FY2=FT2* (BYI+CTI*BYJ) 

FZ2=FT2*  (BZH-CTI*BZJ) 

FX(I)   =FX(I)   -FXl 

FX(I+1)=FX(I+1)+FX1-FX2 

FX ( 1+2 ) =FX ( 1+2 ) +FX2 

FY (I)   =FY(I)   -FYl 

FY (I+l) =FY (I+l) +FY1-FY2 

FY(I+2)=FY(I+2)+FY2 

FZ(I)   =FZ(I)   -FZl 

FZ(I+1)=FZ(I+1)+FZ1-FZ2 

FZ(I+2)=FZ(I+2)+FZ2 
C 

C    CALCULATION  OF  TORSION  FORCES 
C    CHECK  WHETHER  FIRST  ATOM  OF  MOLECULE 

IF(M0D(I+NAM-1,NAM) .EQ.O)  GOTO  600 
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CALL  RYF0R(BX(I-1) ,BY(I-1) ,BZ(I-1) , BX (I) , BY (I) , BZ (I) , 
&    BX(I+1) ,BY(I+1) ,BZ(I+1) ,FX1,FY1,FZ1,FX2,FY2,FZ2, 
&    FX3,FY3,FZ3,FX4,FY4,FZ4,ETQ,CODI,GTHE(I-l) ) 
C 

I COUNT= I COUNT+ 1 

IF(CODI.LT.-0.50)  ITRAN=ITRAN+1 

FX(I-1)=FX(I-1)+FX1 

FY(I-1)=FY(I-1)+FY1 

FZ(I-1)=FZ(I-1)+FZ1 

FX(I)   =FX(I)   +FX2 

FY (I)   =FY(I)   +FY2 

FZ(I)   =FZ(I)   +FZ2 

FX(I+1)=FX(I+1)+FX3 

FY(I+1)=FY(I+1)+FY3 

FZ (I+l) =FZ (I+l) +FZ3 

FX ( 1+2 ) =FX ( 1+2 ) +FX4 

FY ( 1+2 ) =FY ( 1+2 ) +FY4 

FZ(I+2)=FZ(I+2)+FZ4 

ETOR=ETOR+ETQ 
600     CONTINUE 
C 
C   CALCULATE  GROUP-SHELL  INTERACTION 

EWALL-0 . 0 

RCHKSQ= ( RWALL- 2 . 0 ) *  *  2 

DO  650  1=1, NP 
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ICHK=I-H-NAM 
C 

XI=XO(I) 
YI=YO(I) 
ZI=ZO(I) 

c 

C   DISTANCE  FROM  CENTER  TO  GROUP 

RSQ=XI*XI+YI*YI+ZI*ZI 
C   IF  HEAD  GROUP,  CALCULATE  INTERACTIONS  WITH  WALL  USING 
C      A  POTENTIAL  OF  THE  FORM  —   GAMMA* (R-RO) **2 

IF(MOD(ICHK,NAM) .NE.O)  GO  TO  384 
C   DURING  SCALE-DOWN,  DO  NOT  CALCULATE  HEAD-WALL  INTERACTION 
C      GOTO  650 
C 

RR=SQRT(RSQ) 

RI=1.0/RR 

RD=RWALL-RR 
C 

EWALL=EWALL+FORCON* (RD-REQUIL) * (RD-REQUIL) 

FWALL=-2 . 0*FORCON* (RD-REQUIL) 

FTOTWL=FTOTWL+FWALL 
C 

FX ( I ) =FX ( I ) -FWALL*XI *RI 

FY ( I ) =FY ( I ) -FWALL* YI *RI 

FZ (I) =FZ (I) -FWALL*ZI*RI 
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GO  TO  650 
C  IF  NOT  A  HEAD  GROUP,  CALCULATE  REPULSIVE  INTERACTION 
C      WITH  THE  WALL 

384      IF(RSQ.LT.RCHKSQ)  GOTO  650 
C 
C   DISTANCE  FROM  GROUP  TO  SHELL 

RR=SQRT(RSQ) 

RI=1.0/RR 

RD=RWALL-RR 

RDR=1.0/RD 

R3  =RDR*RDR*RDR 

R6  =R3*R3 

R12=R6*R6 
C 
C   FORCE  AND  ENERGY  FOR  R-12  WALL 

EWALL=EWALL+R12 

FWALL=12 . 0*R12*RDR 

FTOTWL=FTOTWL+FWALL 

FX ( I ) =FX ( I ) -FWALL*XI *RI 

FY ( I ) =FY ( I ) -FWALL* YI*RI 

FZ (I) =FZ (I) -FWALL*ZI*RI 
650    CONTINUE 
C 

C   CALCULATE  WALL  PRESSURE  IN  ATMOSPHERES 
EXTPRS=0. 0085382 *FT0TWL/RWALL**2 
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C 

C   CALCULATE  INTRAMOLECULAR  AND  TOTAL  ENERGY 
EINTR=ETOR  +  EBON  +  EBEN 
TOTE  =  TOTLJ  +  EINTR  +  EWALL 
TRFRAC=FLOAT ( ITRAN) /FLOAT ( I COUNT) 

RETURN 

END 


SUBROUTINE  RYFOR ( BONXl , BONYl , BONZ 1 , B0NX2 , B0NY2 , 
&   BONZ 2 , B0NX3 , B0NY3 , BONZ 3 , FORXl , FORYl , FORZl , F0RX2 , 
&   F0RY2 , F0RZ2 , F0RX3 , F0RY3 , F0RZ3 , F0RX4 , F0RY4 , F0RZ4 , 
&   ET0R,C0DI,GR) 
C 

DATA  A1,A2,A3,A4,A5/-1.4  62,-1.578,0.3  68, 
&       3.156,3.788/ 
C 

C   CALCULATE  NORMAL  TO  PLANE  OF  BONDS  1  AND  2 
E=BONYl*BONZ2-BONZl*BONY2 
D=BONZl*BONX2-BONXl*BONZ2 
C=BONXl*BONY2-BONYl*BONX2 
AN0RM2=C*C+D*D+E*E 
ANORM  =SQRT(AN0RM2) 
AN0RMR=1 . /ANORM 
AX=E*ANORMR 
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AY=D*ANORMR 

AZ=C*ANORMR 
C 
C   CALCULATE  NORMAL  TO  PLANE  OF  BONDS  2  AND  3 

H=BONY2*BONZ3-BONZ2*BONY3 

G=B0NZ2  *B0NX3 -B0NX2  *BONZ  3 

F=B0NX2  *B0NY3 -B0NY2  *B0NX3 

BN0RM2=F*F+G*G+H*H 

BNORM  =SQRT(BN0RM2) 

BNORMR= 1 . / BNORM 

BX=H*BNORMR 

BY=G*BNORMR 

BZ=F*BNORMR 
C 
C   DETERMINE  COSINE  OF  THE  DIHEDRAL  ANGLE,  PHI 

CODI=  AX*BX+AY*BY+AZ*BZ 
C 
C   DETERMINE  TORSIONAL  ENERGY 

CODISQ=CODI*CODI 

C0DI3  =CODI*CODISQ 

ETOR  .=  GR*(1.116+Al*CODI+A2*CODISQ+A3*CODI3 
&  +A4*CODISQ*CODISQ+A5*CODISQ*CODI3 ) 

C 
C   CALCULATE  DU/COS(PHI) 

DUDCO= (Al+2 . *A2*CODI+3 . *A3*CODISQ+4 . *A4*C0DISQ* 
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&    CODH-5.*A5*CODISQ*CODISQ)  *GR 

C 

C   START  EVALUATION  OF  FORCE  ON  GROUP  1 

C     CALCULATE  DANORM/DRl 

FACT1= (C*BONY2-D*BONZ2 ) /AN0RM2 
DADX1X=  -FACT1*AX 
DADX1Y=  -FACT1*AY-B0NZ2*AN0RMR 
DADX1Z=  -FACT1*AZ+B0NY2*AN0RMR 

C 

FACT2  =(E*BONZ2-C*BONX2)/ANORM2 
DADY1X=  -FACT2*AX+BONZ2*ANORMR 
DADY1Y=  -FACT2*AY 
DADY1Z=  -FACT2*AZ-BONX2*ANORMR 

C 

FACT3  =(D*BONX2-E*BONY2)/ANORM2 
DADZ1X=  -FACT3*AX-BONY2*ANORMR 
DADZ1Y=  -FACT3*AY+BONX2*ANORMR 
DADZ1Z=  -FACTS *AZ 

C 

C   CALCULATE  D(AB)/DR1 

DABDX1=BX*DADX1X+BY*DADX1Y+BZ*DADX1Z 
DABDY1=BX*DADY1X+BY*DADY1Y+BZ*DADY1Z 
DABDZ1=BX*DADZ1X+BY*DADZ1Y+BZ*DADZ1Z 

C 

C   CALCULATE  FORCE  ON  GROUP  1 
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F0RX1=  -DUDC0*DABDX1 

F0RY1=  -DUDC0*DABDY1 

F0RZ1=  -DUDC0*DABDZ1 
C 

C   START  EVALUATION  OF  FORCE  ON  GROUP  2 
C     CALCULATE  DAN0RM/DR2 

FACT4  =(D*B0NZ1-C*B0NY1)/AN0RM2 

DADX2X=  -FACT4*AX-DADX1X 

DADX2Y=  -FACT4*AY+BONZl*ANORMR-DADXlY 

DADX2Z=  -FACT4*AZ-B0NY1*AN0RMR-DADX1Z 
C 

FACTS  =(C*B0NX1-E*B0NZ1)/AN0RM2 

DADY2X=  -FACT5*AX-B0NZ1*AN0RMR-DADY1X 

DADY2Y=  -FACT5*AY-DADY1Y 

DADY2Z=  -FACT5*AZ+B0NX1*AN0RMR-DADY1Z 
C 

FACT6  =(E*B0NY1-D*B0NX1)/AN0RM2 

DADZ2X=  -FACT6*AX+B0NY1*AN0RMR-DADZ1X 

DADZ2Y=  -FACT6*AY-B0NX1*AN0RMR-DADZ1Y 

DADZ2Z=  -FACT6*AZ-DADZ1Z 
C 
C   EVALUATE  D(BN0RM)/DR2 

FACT7  =(F*BONY3-G*BONZ3)/BNORM2 

DBDX2X=  -FACT7*BX 

DBDX2Y=  -FACT7*BY-BONZ3*BNORMR 
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DBDX2Z=  -FACT7*BZ+BONY3*BNORMR 
C 

FACTS  =(H*BONZ3-F*BONX3)/BNORM2 

DBDY2X=  -FACT8*BX+BONZ3*BNORMR 

DBDY2Y=  -FACT8*BY 

DBDY2Z=  -FACT8*BZ-BONX3*BNORMR 
C 

FACT9  =(G*BONX3-H*BONY3)/BNORM2 

DBDZ2X=  -FACT9*BX-BONY3*BNORMR 

DBDZ2Y=  -FACT9*BY+BONX3*BNORMR 

DBDZ2Z=  -FACT9*BZ 
C 
C   CALCULATE  D  C0S(PHI)/DR2 

DABDX2=BX*DADX2X+BY*DADX2Y+BZ*DADX2Z 
&  +AX*DBDX2X+AY*DBDX2Y+AZ*DBDX2Z 

DABDY2=BX*DADY2X+BY*DADY2Y+BZ*DADY2Z 
&         +AX*DBDY2X+AY*DBDY2Y+AZ*DBDY2Z 

DABDZ2=BX*DADZ2X+BY*DADZ2Y+BZ*DADZ2Z 
&         +AX*DBDZ2X+AY*DBDZ2Y+AZ*DBDZ2Z 
C 
C   FORCE  ON  GROUP  2 

F0RX2  =  -DUDC0*DABDX2 

F0RY2  =  -DUDC0*DABDY2 

F0RZ2  =  -DUDC0*DABDZ2 
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C   START  EVALUATION  OF  FORCE  ON  GROUP  4 

C     CALCULATE  D(BN0RM)/DR4 

FACTl  =(F*BONY2-G*BONZ2)/BNORM2 
DBDX4X=  -FACT1*BX 
DBDX4Y=  -FACT1*BY-B0NZ2*BN0RMR 
DBDX4Z=  -FACT1*BZ+B0NY2*BN0RMR 

C 

FACT2  =(H*BONZ2-F*BONX2)/BNORM2 
DBDY4X=  -FACT2*BX+BONZ2*BNORMR 
DBDY4Y=  -FACT2*BY 
DBDY4Z=  -FACT2*BZ-BONX2*BNORMR 

C 

FACT3  =(G*BONX2-H*BONY2)/BNORM2 
DBDZ4X=  -FACT3*BX-BONY2*BNORMR 
DBDZ4Y=  -FACT3*BY+BONX2*BNORMR 
DBDZ4Z=  -FACT3*BZ 

C 

C   EVALUATE  D  C0S(PHI)/DR4 

DABDX4=AX*DBDX4X+AY*DBDX4Y+AZ*DBDX4Z 
DABDY4=AX*DBDY4X+AY*DBDY4Y+AZ*DBDY4Z 
DABDZ4=AX*DBDZ4X+AY*DBDZ4Y+AZ*DBDZ4Z 

C 

C  EVALUATE  FORCE  ON  GROUP  4 
F0RX4=  -DUDC0*DABDX4 
F0RY4=  -DUDC0*DABDY4 
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F0RZ4=  -DUDC0*DABDZ4 
C 
C   EVALUATE  FORCE  ON  GROUP  3  BY  NEWTON'S  THIRD  LAW 

F0RX3  =  -FORX1-FORX2-FORX4 

F0RY3  =  -FORY1-FORY2-FORY4 

F0RZ3  =  -FORZ1-FORZ2-FORZ4 
C 

RETURN 

END 


BLOCK  DATA 

COMMON/DER/X2(216) , Y2 (216) , Z2 (216) , X3 (216) , Y3 (216) , 
&   Z3(216) ,X4(216) ,Y4(216) ,Z4(216) ,X5(216) ,Y5(216) , 
&   Z5(216) 

COMMON/FOR/FX(216) , FY (216) , FZ (216) 

COMMON/DISP/DAX(216) , DAY (216) ,DAZ(216) ,X0L(216) , 
&   Y0L(216) ,Z0L(216) 

COMMON/ PRO P/SUME , XSUM 

COMMON/MASS/MASS (216) 

COMMON/ INTRA/BO (216) ,GAB(216) ,CT0(216) ,GTHA(216) , 
&   GTHE(216) 

DATA  X3,Y3,Z3,X4,Y4,Z4,X5,Y5,Z5/1944*0./ 

DATA  DAX,DAY,DAZ,FX,FY,FZ/1296*0./ 

DATA  MASS/216*1/ 
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C    FOR  SCALE-DOWN,  GTHE  SET  TO  GTHE/10. 

DATA  B0,GAB,CT0,GTHA,GTHE/216*.3  8475,216*3  53  22.1957, 
&   216*-. 37703  2  668,216*310.26253,216*19.8424  8/ 
C    NORMAL  GTHE  IS  19.84248 

END 
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